Changeset 11639 for branches/UKMO

Timestamp:

2019-10-02T13:13:03+02:00 (5 years ago)

Author:

rrenshaw

Message:

regional mean and transport diagnostics added for reanalysis

Location:

branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO

Files:

• OPA_SRC/DIA/diaar5.F90 (modified) (10 diffs)
• OPA_SRC/DIA/diadct.F90 (modified) (41 diffs)
• OPA_SRC/DIA/diapea.F90 (added)
• OPA_SRC/DIA/diaregmean.F90 (added)
• OPA_SRC/DIA/diatmb.F90 (modified) (1 diff)
• OPA_SRC/DIA/diawri.F90 (modified) (2 diffs)
• OPA_SRC/IOM/iom.F90 (modified) (3 diffs)
• OPA_SRC/OBS/obs_level_search.h90 (modified) (2 diffs)
• OPA_SRC/OBS/obsinter_z1d.h90 (modified) (3 diffs)
• OPA_SRC/nemogcm.F90 (modified) (3 diffs)
• TOP_SRC/trcbdy.F90 (modified) (1 diff)

branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90

@@ -39 +39 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:  ) ::   thick0       ! ocean thickness (interior domain)
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sn0          ! initial salinity
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tn0          ! initial temperature
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshthster_mat         ! ssh_thermosteric height
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshhlster_mat         ! ssh_halosteric height
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sshsteric_mat         ! ssh_steric height
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   zbotpres_mat          ! bottom pressure
       
    !! * Substitutions
@@ -56 +61 @@
       !!----------------------------------------------------------------------
       !
-      ALLOCATE( area(jpi,jpj), thick0(jpi,jpj) , sn0(jpi,jpj,jpk) , STAT=dia_ar5_alloc )
+      ALLOCATE( area(jpi,jpj), thick0(jpi,jpj) , sn0(jpi,jpj,jpk), tn0(jpi,jpj,jpk) , &
+          & sshthster_mat(jpi,jpj),sshhlster_mat(jpi,jpj),sshsteric_mat(jpi,jpj), &
+          & zbotpres_mat(jpi,jpj),STAT=dia_ar5_alloc )
       !
       IF( lk_mpp             )   CALL mpp_sum ( dia_ar5_alloc )
@@ -85 +92 @@
       CALL wrk_alloc( jpi , jpj , jpk        , zrhd      , zrhop    )
       CALL wrk_alloc( jpi , jpj , jpk , jpts , ztsn                 )
+      
+      sshthster_mat(:,:) = 0._wp  
+      sshhlster_mat(:,:) = 0._wp  
+      sshsteric_mat(:,:) = 0._wp  
+      zbotpres_mat(:,:)  = 0._wp  
 
       zarea_ssh(:,:) = area(:,:) * sshn(:,:)
@@ -121 +133 @@
       zssh_steric = - zarho / area_tot
       CALL iom_put( 'sshthster', zssh_steric )
+      sshthster_mat(:,:) =  -zbotpres(:,:)
+      CALL iom_put( 'sshthster_mat', sshthster_mat )
+
+      !                     
+      ztsn(:,:,:,jp_tem) = tn0(:,:,:)                    ! thermohaline ssh
+      ztsn(:,:,:,jp_sal) = tsn(:,:,:,jp_sal) 
+      CALL eos( ztsn, zrhd, fsdept_n(:,:,:) )                       ! now in situ density using initial temperature
+      !
+      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
+            DO ji=1,jpi
+               DO jj=1,jpj
+                  zbotpres(ji,jj) = zbotpres(ji,jj) + sshn(ji,jj) * zrhd(ji,jj,mikt(ji,jj)) + riceload(ji,jj)
+               END DO
+            END DO
+         ELSE
+            zbotpres(:,:) = zbotpres(:,:) + sshn(:,:) * zrhd(:,:,1)
+         END IF
+      END IF
+      !                                         
+      zarho = SUM( area(:,:) * zbotpres(:,:) ) 
+      IF( lk_mpp )   CALL mpp_sum( zarho )
+      zssh_steric = - zarho / area_tot
+      CALL iom_put( 'sshhlster', zssh_steric )
+      sshhlster_mat(:,:) = -zbotpres(:,:)
+      CALL iom_put( 'sshhlster_mat', sshhlster_mat )
+      
+
       
       !                                         ! steric sea surface height
@@ -147 +191 @@
       zssh_steric = - zarho / area_tot
       CALL iom_put( 'sshsteric', zssh_steric )
+      sshsteric_mat(:,:) = -zbotpres(:,:) 
+      CALL iom_put( 'sshsteric_mat', sshsteric_mat )
       
       !                                         ! ocean bottom pressure
       zztmp = rau0 * grav * 1.e-4_wp               ! recover pressure from pressure anomaly and cover to dbar = 1.e4 Pa
       zbotpres(:,:) = zztmp * ( zbotpres(:,:) + sshn(:,:) + thick0(:,:) )
-      CALL iom_put( 'botpres', zbotpres )
-
+      zbotpres_mat(:,:) = zbotpres(:,:)
+      CALL iom_put( 'botpres', zbotpres_mat )
+      
       !                                         ! Mean density anomalie, temperature and salinity
       ztemp = 0._wp
@@ -211 +258 @@
       REAL(wp) ::   zztmp  
       REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   zsaldta   ! Jan/Dec levitus salinity
+      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   ztemdta   ! Jan/Dec levitus salinity
       ! reading initial file
       LOGICAL  ::   ln_tsd_init      !: T & S data flag
@@ -234 +282 @@
       !
       CALL wrk_alloc( jpi , jpj , jpk, jpts, zsaldta )
+      CALL wrk_alloc( jpi , jpj , jpk, jpts, ztemdta )
       !                                      ! allocate dia_ar5 arrays
       IF( dia_ar5_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'dia_ar5_init : unable to allocate arrays' )
@@ -253 +302 @@
       CALL iom_get  ( inum, jpdom_data, TRIM(sn_sal%clvar), zsaldta(:,:,:,2), 12 )
       CALL iom_close( inum )
+
+      CALL iom_open ( TRIM( cn_dir )//TRIM(sn_tem%clname), inum )
+      CALL iom_get  ( inum, jpdom_data, TRIM(sn_tem%clvar), ztemdta(:,:,:,1), 1  )
+      CALL iom_get  ( inum, jpdom_data, TRIM(sn_tem%clvar), ztemdta(:,:,:,2), 12 )
+      CALL iom_close( inum )
       sn0(:,:,:) = 0.5_wp * ( zsaldta(:,:,:,1) + zsaldta(:,:,:,2) )        
-      sn0(:,:,:) = sn0(:,:,:) * tmask(:,:,:)
+      tn0(:,:,:) = 0.5_wp * ( ztemdta(:,:,:,1) + ztemdta(:,:,:,2) )        
+      sn0(:,:,:) = sn0(:,:,:) * tmask(:,:,:)  
+      tn0(:,:,:) = tn0(:,:,:) * tmask(:,:,:)
       IF( ln_zps ) THEN               ! z-coord. partial steps
          DO jj = 1, jpj               ! interpolation of salinity at the last ocean level (i.e. the partial step)
@@ -262 +318 @@
                   zztmp = ( gdept_1d(ik) - gdept_0(ji,jj,ik) ) / ( gdept_1d(ik) - gdept_1d(ik-1) )
                   sn0(ji,jj,ik) = ( 1._wp - zztmp ) * sn0(ji,jj,ik) + zztmp * sn0(ji,jj,ik-1)
+                  tn0(ji,jj,ik) = ( 1._wp - zztmp ) * tn0(ji,jj,ik) + zztmp * tn0(ji,jj,ik-1)
                ENDIF
             END DO
@@ -268 +325 @@
       !
       CALL wrk_dealloc( jpi , jpj , jpk, jpts, zsaldta )
+      CALL wrk_dealloc( jpi , jpj , jpk, jpts, ztemdta )
       !
       IF( nn_timing == 1 )   CALL timing_stop('dia_ar5_init')

branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90

@@ -35 +35 @@
   USE phycst          ! physical constants
   USE in_out_manager  ! I/O manager
+  USE iom             ! I/0 library
   USE daymod          ! calendar
   USE dianam          ! build name of file
@@ -65 +66 @@
   !! * Shared module variables
   LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .TRUE.   !: model-data diagnostics flag
+  LOGICAL, PUBLIC ::   ln_dct_calc_noos_25h !: Calcuate noos 25 h means
+  LOGICAL, PUBLIC ::   ln_dct_calc_noos_hr !: Calcuate noos hourly means
+  ! JT
+  LOGICAL, PUBLIC ::   ln_dct_iom_cont !: Use IOM Output?
+  LOGICAL, PUBLIC ::   ln_dct_ascii    !: Output ascii or binary
+  LOGICAL, PUBLIC ::   ln_dct_h        !: Output hourly instantaneous or mean values
+  ! JT
 
   !! * Module variables
-  INTEGER :: nn_dct      = 1     ! Frequency of computation
-  INTEGER :: nn_dctwri   = 1     ! Frequency of output
-  INTEGER :: nn_secdebug = 0     ! Number of the section to debug
-  INTEGER :: nn_dct_h    = 1     ! Frequency of computation for NOOS hourly files
-  INTEGER :: nn_dctwri_h = 1     ! Frequency of output for NOOS hourly files
+  INTEGER :: nn_dct        ! Frequency of computation
+  INTEGER :: nn_dctwri     ! Frequency of output
+  INTEGER :: nn_secdebug   ! Number of the section to debug
+  INTEGER :: nn_dct_h      ! Frequency of computation for NOOS hourly files
+  INTEGER :: nn_dctwri_h   ! Frequency of output for NOOS hourly files
    
   INTEGER, PARAMETER :: nb_class_max  = 12   ! maximum number of classes, i.e. depth levels or density classes
-  INTEGER, PARAMETER :: nb_sec_max    = 30   ! maximum number of sections
-  INTEGER, PARAMETER :: nb_point_max  = 300  ! maximum number of points in a single section
+  INTEGER, PARAMETER :: nb_sec_max    = 100  ! maximum number of sections
+  INTEGER, PARAMETER :: nb_point_max  = 375  ! maximum number of points in a single section
   INTEGER, PARAMETER :: nb_type_class       = 14   ! types of calculations, i.e. pos transport, neg transport, heat transport, salt transport
   INTEGER, PARAMETER :: nb_3d_vars    = 5
@@ -134 +142 @@
      ALLOCATE(transports_3d_h(nb_3d_vars,nb_sec_max,nb_point_max,jpk), STAT=ierr(3) )
      ALLOCATE(transports_2d_h(nb_2d_vars,nb_sec_max,nb_point_max)    , STAT=ierr(4) )
-     ALLOCATE(z_hr_output(nb_sec_max,168,nb_class_max)               , STAT=ierr(5) ) ! 168 = 24 * 7days
+     ALLOCATE(z_hr_output(nb_sec_max,3,nb_class_max)                , STAT=ierr(5) )
 
      diadct_alloc = MAXVAL( ierr ) 
@@ -149 +157 @@
      !!
      !!---------------------------------------------------------------------
-     NAMELIST/namdct/nn_dct,nn_dctwri,nn_secdebug
-     INTEGER  ::   ios                 ! Local integer output status for namelist read
+     NAMELIST/namdct/nn_dct,ln_dct_h,nn_dctwri,ln_dct_ascii,nn_secdebug,ln_dct_calc_noos_25h,ln_dct_calc_noos_hr,ln_dct_iom_cont
+     INTEGER           ::   ios,jsec        ! Local integer output status for namelist read
+     CHARACTER(len=3)  ::   jsec_str        ! name of the jsec
+     LOGICAL       ::   verbose     
+     verbose = .false.
 
      IF( nn_timing == 1 )   CALL timing_start('dia_dct_init')
@@ -164 +175 @@
 
      IF( ln_NOOS ) THEN
-        nn_dct=3600./rdt         ! hard coded for NOOS transects, to give 25 hour means 
+
+        !Do calculation for daily, 25hourly mean every hour
+        nn_dct=3600./rdt         ! hard coded for NOOS transects, to give 25 hour means from hourly instantaneous values
+
+        !write out daily, 25hourly mean every day
         nn_dctwri=86400./rdt
-        nn_dct_h=1       ! hard coded for NOOS transects, to give hourly data
+
+
+        !nn_dct_h=1       ! hard coded for NOOS transects, to give hourly data    
+        ! If you want hourly instantaneous values, you only do the calculation every 12 timesteps (if rdt = 300)
+        ! and output it every 12 time steps. For this, you set the ln_dct_h to be True, and it calcuates it automatically
+        ! if you want hourly mean values, set ln_dct_h to be False, and it will do the calculate every time step.
+        !
+        !SELECT CASE( ln_dct_h )
+        !   CASE(.TRUE.)  
+        !      nn_dct_h=3600./rdt
+        !   CASE(.FALSE.)
+        !      nn_dct_h=1
+        !END SELECT
+        
+        IF ( ln_dct_h ) THEN
+            nn_dct_h=3600./rdt
+        ELSE
+            nn_dct_h=1.
+        ENDIF    
+        
+        !JT write out hourly calculation every hour
         nn_dctwri_h=3600./rdt
      ENDIF
 
      IF( lwp ) THEN
+       IF( verbose ) THEN
         WRITE(numout,*) " "
         WRITE(numout,*) "diadct_init: compute transports through sections "
         WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
         IF( ln_NOOS ) THEN
+           WRITE(numout,*) "       Calculate NOOS hourly output: ln_dct_calc_noos_hr = ",ln_dct_calc_noos_hr
+           WRITE(numout,*) "       Calculate NOOS 25 hour mean output: ln_dct_calc_noos_hr = ",ln_dct_calc_noos_25h
+           WRITE(numout,*) "       Use IOM Output: ln_dct_iom_cont = ",ln_dct_iom_cont
+           WRITE(numout,*) "       Output in ASCII (True) or Binary (False): ln_dct_ascii = ",ln_dct_ascii
+           WRITE(numout,*) "       Frequency of hourly computation - instantaneous (TRUE) or hourly mean (FALSE): ln_dct_h  = ",ln_dct_h
+
            WRITE(numout,*) "       Frequency of computation hard coded to be every hour: nn_dct    = ",nn_dct
            WRITE(numout,*) "       Frequency of write hard coded to average 25 instantaneous hour values: nn_dctwri = ",nn_dctwri
-           WRITE(numout,*) "       Frequency of hourly computation hard coded to be every timestep: nn_dct_h  = ",nn_dct_h
+           WRITE(numout,*) "       Frequency of hourly computation (timestep) : nn_dct_h  = ",nn_dct_h
+           WRITE(numout,*) "       Frequency of hourly computation Not hard coded to be every timestep, or : nn_dct_h  = ",nn_dct_h
            WRITE(numout,*) "       Frequency of hourly write hard coded to every hour: nn_dctwri_h = ",nn_dctwri_h
         ELSE
@@ -183 +226 @@
            WRITE(numout,*) "       Frequency of write:       nn_dctwri = ",nn_dctwri
         ENDIF
+       ENDIF
 
         IF      ( nn_secdebug .GE. 1 .AND. nn_secdebug .LE. nb_sec_max )THEN
@@ -195 +239 @@
 
      ENDIF
-
-     !Read section_ijglobal.diadct
-     CALL readsec
+     
+     
+     IF ( ln_NOOS ) THEN
+        IF ( ln_dct_calc_noos_25h .or. ln_dct_calc_noos_hr ) CALL readsec
+     ENDIF
 
      !open output file
-     IF( lwm ) THEN
+     IF( lwp ) THEN
         IF( ln_NOOS ) THEN
-           CALL ctl_opn( numdct_NOOS  ,'NOOS_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
-           CALL ctl_opn( numdct_NOOS_h,'NOOS_transport_h', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
+           WRITE(numout,*) "diadct_init: Open output files. ASCII? ",ln_dct_ascii
+           WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
+           IF  ( ln_dct_ascii ) THEN
+               if ( ln_dct_calc_noos_25h ) CALL ctl_opn( numdct_NOOS  ,'NOOS_transport'  , 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .TRUE. )
+               if ( ln_dct_calc_noos_hr )  CALL ctl_opn( numdct_NOOS_h,'NOOS_transport_h', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .TRUE. )
+           ELSE
+               if ( ln_dct_calc_noos_25h ) CALL ctl_opn( numdct_NOOS  ,'NOOS_transport_bin'  , 'REPLACE', 'UNFORMATTED', 'SEQUENTIAL', -1, numout,  .TRUE. )
+               if ( ln_dct_calc_noos_hr )  CALL ctl_opn( numdct_NOOS_h,'NOOS_transport_bin_h', 'REPLACE', 'UNFORMATTED', 'SEQUENTIAL', -1, numout,  .TRUE. )
+           ENDIF
         ELSE
-           CALL ctl_opn( numdct_vol,  'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
-           CALL ctl_opn( numdct_heat, 'heat_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
-           CALL ctl_opn( numdct_salt, 'salt_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
+           CALL ctl_opn( numdct_vol,  'volume_transport', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
+           CALL ctl_opn( numdct_heat, 'heat_transport'  , 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
+           CALL ctl_opn( numdct_salt, 'salt_transport'  , 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
         ENDIF
      ENDIF
 
      ! Initialise arrays to zero 
-     transports_3d(:,:,:,:)  =0.0 
-     transports_2d(:,:,:)    =0.0 
-     transports_3d_h(:,:,:,:)=0._wp
-     transports_2d_h(:,:,:)  =0._wp
-     z_hr_output(:,:,:)      =0._wp
+     transports_3d(:,:,:,:)   =0._wp
+     transports_2d(:,:,:)     =0._wp
+     transports_3d_h(:,:,:,:) =0._wp
+     transports_2d_h(:,:,:)   =0._wp
+     z_hr_output(:,:,:)       =0._wp
 
      IF( nn_timing == 1 )   CALL timing_stop('dia_dct_init')
+
+     IF (ln_dct_iom_cont) THEN
+        IF( lwp ) THEN
+            WRITE(numout,*) " "
+            WRITE(numout,*) "diadct_init: using xios iom_put for output: field_def.xml and iodef.xml code"
+            WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
+            WRITE(numout,*) ""
+            WRITE(numout,*) "      field_def.xml"
+            WRITE(numout,*) "      ~~~~~~~~~~~~~"
+            WRITE(numout,*) ""
+            WRITE(numout,*) ""
+      
+            WRITE(numout,*)  '      <field_group id="noos_cross_section" domain_ref="1point" axis_ref="noos" operation="average">'
+
+            DO jsec=1,nb_sec
+               WRITE (jsec_str, "(I3.3)") jsec
+               
+               WRITE(numout,*)  '          <field id="noos_'//jsec_str//'_trans" long_name="' // TRIM(secs(jsec)%name) // ' 25h mean NOOS transport cross-section number: '//jsec_str//' (total, positive, negative)" unit="m^3/s" />'
+               WRITE(numout,*)  '          <field id="noos_'//jsec_str//'_heat" long_name="' // TRIM(secs(jsec)%name) // ' 25h mean NOOS heat cross-section number: '//jsec_str//' (total, positive, negative)" unit="J/s" />'
+               WRITE(numout,*)  '          <field id="noos_'//jsec_str//'_salt" long_name="' // TRIM(secs(jsec)%name) // ' 25h mean NOOS salt cross-section number: '//jsec_str//' (total, positive, negative)" unit="g/s" />'
+               
+            ENDDO
+            
+            WRITE(numout,*)  '      </field_group>'
+            
+            WRITE(numout,*) ""
+            WRITE(numout,*) ""
+            WRITE(numout,*) "      iodef.xml"
+            WRITE(numout,*) "      ~~~~~~~~~"
+            WRITE(numout,*) ""
+            WRITE(numout,*) ""
+            
+            WRITE(numout,*)  '      <file_group id="1d" output_freq="1d" output_level="10" enabled=".TRUE.">'
+            WRITE(numout,*) ""
+            WRITE(numout,*)  '          <file id="noos_cross_section" name="NOOS_transport">'
+            DO jsec=1,nb_sec
+               WRITE (jsec_str, "(I3.3)") jsec
+               
+               WRITE(numout,*)  '              <field field_ref="noos_'//jsec_str//'_trans" />'
+               WRITE(numout,*)  '              <field field_ref="noos_'//jsec_str//'_heat" />'
+               WRITE(numout,*)  '              <field field_ref="noos_'//jsec_str//'_salt" />'
+               
+            ENDDO
+            WRITE(numout,*)  '          </file>'
+            WRITE(numout,*) ""
+            WRITE(numout,*)  '      </file_group>'
+            
+            WRITE(numout,*) ""
+            WRITE(numout,*) ""
+            WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
+            WRITE(numout,*) ""
+        
+        ENDIF
+     ENDIF
+
+
      !
   END SUBROUTINE dia_dct_init
@@ -231 +340 @@
      !!  Method  :: All arrays initialised to zero in dct_init 
      !!             Each nn_dct time step call subroutine 'transports' for 
-     !!               each section to sum the transports over each grid cell. 
+     !!               each section to sum the transports.
      !!             Each nn_dctwri time step: 
      !!               Divide the arrays by the number of summations to gain 
@@ -252 +361 @@
      REAL(wp), POINTER, DIMENSION(:)    :: zwork !   "  
      REAL(wp), POINTER, DIMENSION(:,:,:):: zsum  !   "
+     LOGICAL       ::   verbose     
+     verbose = .false.
 
      !!---------------------------------------------------------------------    
@@ -266 +377 @@
      zsum(:,:,:) = 0.0
 
-     IF( lwp .AND. kt==nit000+nn_dct-1 ) THEN
+     IF( lwp .AND. kt==nit000+nn_dct-1 .AND. verbose ) THEN
          WRITE(numout,*) " "
          WRITE(numout,*) "diadct: compute transport"
          WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~~~~~"
          WRITE(numout,*) "nb_sec = ",nb_sec
+         WRITE(numout,*) "nn_dct = ",nn_dct
+         WRITE(numout,*) "ln_NOOS = ",ln_NOOS
+         WRITE(numout,*) "nb_sec = ",nb_sec
+         WRITE(numout,*) "nb_sec_max = ",nb_sec_max
+         WRITE(numout,*) "nb_type_class = ",nb_type_class
+         WRITE(numout,*) "nb_class_max = ",nb_class_max
      ENDIF
 
- 
-     ! Compute transport and write only at nn_dctwri
-     IF( MOD(kt,nn_dct)==0 .or. &                ! compute transport every nn_dct time steps
-         (ln_NOOS .and. kt==nn_it000 ) )  THEN   ! also include first time step when calculating NOOS 25 hour averages
-
-        DO jsec=1,nb_sec
-
-           !debug this section computing ?
-           IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND.  kt==nit000+nn_dct-1 .AND. lwp ) lldebug=.TRUE. 
-
-           !Compute transport through section  
-           CALL transport(secs(jsec),lldebug,jsec) 
-
-        ENDDO
-             
-        IF( MOD(kt,nn_dctwri)==0 )THEN
-
-           IF( lwp .AND. kt==nit000+nn_dctwri-1 )WRITE(numout,*)"      diadct: average transports and write at kt = ",kt         
-  
-           !! divide arrays by nn_dctwri/nn_dct to obtain average 
-           transports_3d(:,:,:,:)=transports_3d(:,:,:,:)/(nn_dctwri/nn_dct) 
-           transports_2d(:,:,:)  =transports_2d(:,:,:)  /(nn_dctwri/nn_dct) 
- 
-           ! Sum over each class 
-           DO jsec=1,nb_sec 
-              CALL dia_dct_sum(secs(jsec),jsec) 
-           ENDDO 
-
-           !Sum on all procs 
-           IF( lk_mpp )THEN
-              zsum(:,:,:)=0.0_wp
-              ish(1)  =  nb_sec_max*nb_type_class*nb_class_max 
-              ish2    = (/nb_sec_max,nb_type_class,nb_class_max/)
-              DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport(:,:) ; ENDDO
-              zwork(:)= RESHAPE(zsum(:,:,:), ish )
-              CALL mpp_sum(zwork, ish(1))
-              zsum(:,:,:)= RESHAPE(zwork,ish2)
-              DO jsec=1,nb_sec ; secs(jsec)%transport(:,:) = zsum(jsec,:,:) ; ENDDO
-           ENDIF
-
-           !Write the transport
-           DO jsec=1,nb_sec
-
-              IF( lwm .and. .not. ln_NOOS )CALL dia_dct_wri(kt,jsec,secs(jsec))
-              IF( lwm .and.       ln_NOOS )CALL dia_dct_wri_NOOS(kt,jsec,secs(jsec))   ! use NOOS specific formatting
+    
+    IF ( ln_dct_calc_noos_25h ) THEN
+ 
+        ! Compute transport and write only at nn_dctwri
+        IF ( MOD(kt,nn_dct)==0 .or. &               ! compute transport every nn_dct time steps
+            (ln_NOOS .and. kt==nn_it000 ) )  THEN   ! also include first time step when calculating NOOS 25 hour averages
             
-              !nullify transports values after writing
-              transports_3d(:,jsec,:,:)=0.
-              transports_2d(:,jsec,:  )=0.
-              secs(jsec)%transport(:,:)=0. 
-              IF ( ln_NOOS ) CALL transport(secs(jsec),lldebug,jsec)  ! reinitialise for next 25 hour instantaneous average (overlapping values)
-
-           ENDDO
-
-        ENDIF 
-
+            
+
+            DO jsec=1,nb_sec
+
+              lldebug=.FALSE.
+              IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND.  kt==nit000+nn_dct-1 .AND. lwp ) lldebug=.TRUE. 
+
+              !Compute transport through section  
+              CALL transport(secs(jsec),lldebug,jsec) 
+              
+
+            ENDDO
+                
+            IF( MOD(kt,nn_dctwri)==0 )THEN
+            
+            
+
+              IF( lwp .AND. kt==nit000+nn_dctwri-1 .AND. verbose ) WRITE(numout,*)"      diadct: average and write at kt = ",kt
+
+
+              ! Not 24 values, but 25! divide by ((nn_dctwri/nn_dct) +1)
+              !! divide arrays by nn_dctwri/nn_dct  to obtain average
+              transports_3d(:,:,:,:)= transports_3d(:,:,:,:)/((nn_dctwri/nn_dct)+1.)
+              transports_2d(:,:,:)  = transports_2d(:,:,:)  /((nn_dctwri/nn_dct)+1.)
+
+              ! Sum over each class
+              DO jsec=1,nb_sec
+                  CALL dia_dct_sum(secs(jsec),jsec)
+              ENDDO
+    
+              !Sum on all procs 
+              IF( lk_mpp )THEN
+                  zsum(:,:,:)=0.0_wp
+                  ish(1)  =  nb_sec_max*nb_type_class*nb_class_max 
+                  ish2    = (/nb_sec_max,nb_type_class,nb_class_max/)
+                  DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport(:,:) ; ENDDO
+                  zwork(:)= RESHAPE(zsum(:,:,:), ish )
+                  CALL mpp_sum(zwork, ish(1))
+                  zsum(:,:,:)= RESHAPE(zwork,ish2)
+                  DO jsec=1,nb_sec ; secs(jsec)%transport(:,:) = zsum(jsec,:,:) ; ENDDO
+              ENDIF
+
+              !Write the transport
+              DO jsec=1,nb_sec
+
+                  IF( lwp .and. .not. ln_NOOS )CALL dia_dct_wri(kt,jsec,secs(jsec))
+                  !IF( lwp .and.       ln_NOOS )CALL dia_dct_wri_NOOS(kt,jsec,secs(jsec))   ! use NOOS specific formatting
+                  IF(  ln_NOOS )CALL dia_dct_wri_NOOS(kt,jsec,secs(jsec))   ! use NOOS specific formatting
+                
+                
+                  !nullify transports values after writing
+                  transports_3d(:,jsec,:,:)=0.0
+                  transports_2d(:,jsec,:  )=0.0
+                  secs(jsec)%transport(:,:)=0.  
+                  
+                  
+                  IF ( ln_NOOS ) CALL transport(secs(jsec),lldebug,jsec)  ! reinitialise for next 25 hour instantaneous average (overlapping values)
+
+
+
+              ENDDO
+
+            ENDIF 
+
+        ENDIF
+        
+    ENDIF
+    IF ( ln_dct_calc_noos_hr ) THEN
+        IF ( MOD(kt,nn_dct_h)==0 ) THEN            ! compute transport every nn_dct_h time steps
+
+            DO jsec=1,nb_sec
+
+              lldebug=.FALSE.
+              IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND.  kt==nit000+nn_dct_h-1 .AND. lwp ) lldebug=.TRUE. 
+
+              !Compute transport through section  
+              CALL transport_h(secs(jsec),lldebug,jsec) 
+
+            ENDDO
+                
+            IF( MOD(kt,nn_dctwri_h)==0 )THEN
+
+              IF( lwp .AND. kt==nit000+nn_dctwri_h-1 .AND. verbose )WRITE(numout,*)"      diadct: average and write hourly files at kt = ",kt         
+      
+              !! divide arrays by nn_dctwri/nn_dct to obtain average
+                !
+                ! JT - I think this is wrong. I think it is trying to sum over 25 hours, but only dividing by 24.
+                ! I think it might work for daily cycles, but not for monthly cycles,
+                !
+              transports_3d_h(:,:,:,:)=transports_3d_h(:,:,:,:)/(nn_dctwri_h/nn_dct_h)
+              transports_2d_h(:,:,:)  =transports_2d_h(:,:,:)  /(nn_dctwri_h/nn_dct_h)
+
+              ! Sum over each class
+              DO jsec=1,nb_sec
+                  CALL dia_dct_sum_h(secs(jsec),jsec)
+              ENDDO
+
+              !Sum on all procs 
+              IF( lk_mpp )THEN
+                  ish(1)  =  nb_sec_max*nb_type_class*nb_class_max 
+                  ish2    = (/nb_sec_max,nb_type_class,nb_class_max/)
+                  DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport_h(:,:) ; ENDDO
+                  zwork(:)= RESHAPE(zsum(:,:,:), ish )
+                  CALL mpp_sum(zwork, ish(1))
+                  zsum(:,:,:)= RESHAPE(zwork,ish2)
+                  DO jsec=1,nb_sec ; secs(jsec)%transport_h(:,:) = zsum(jsec,:,:) ; ENDDO
+              ENDIF
+
+              !Write the transport
+              DO jsec=1,nb_sec
+
+                  IF( lwp .and. ln_NOOS ) THEN
+                    CALL dia_dct_wri_NOOS_h(kt/nn_dctwri_h,jsec,secs(jsec))   ! use NOOS specific formatting
+                  endif
+                  !nullify transports values after writing
+                  transports_3d_h(:,jsec,:,:)=0.0
+                  transports_2d_h(:,jsec,:)=0.0
+                  secs(jsec)%transport_h(:,:)=0.0
+                  
+                  ! for hourly mean or hourly instantaneous, you don't initialise! start with zero!
+                  !IF ( ln_NOOS ) CALL transport_h(secs(jsec),lldebug,jsec)  ! reinitialise for next 25 hour instantaneous average (overlapping values)
+
+              ENDDO
+
+            ENDIF 
+
+        ENDIF    
+        
      ENDIF
-
-     IF ( MOD(kt,nn_dct_h)==0 ) THEN            ! compute transport every nn_dct_h time steps
-
-        DO jsec=1,nb_sec
-
-           !lldebug=.FALSE.
-           IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND.  kt==nit000+nn_dct_h-1 .AND. lwp ) lldebug=.TRUE. 
-
-           !Compute transport through section  
-           CALL transport_h(secs(jsec),lldebug,jsec) 
-
-        ENDDO
-             
-        IF( MOD(kt,nn_dctwri_h)==0 )THEN
-
-           IF( lwp .AND. kt==nit000+nn_dctwri_h-1 )WRITE(numout,*)"      diadct: average and write hourly files at kt = ",kt         
-  
-           !! divide arrays by nn_dctwri/nn_dct to obtain average
-           transports_3d_h(:,:,:,:)=transports_3d_h(:,:,:,:)/(nn_dctwri_h/nn_dct_h)
-           transports_2d_h(:,:,:)  =transports_2d_h(:,:,:)  /(nn_dctwri_h/nn_dct_h)
-
-           ! Sum over each class
-           DO jsec=1,nb_sec
-              CALL dia_dct_sum_h(secs(jsec),jsec)
-           ENDDO
- 
-           !Sum on all procs 
-          IF( lk_mpp )THEN
-              ish(1)  =  nb_sec_max*nb_type_class*nb_class_max 
-              ish2    = (/nb_sec_max,nb_type_class,nb_class_max/)
-              DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport_h(:,:) ; ENDDO
-              zwork(:)= RESHAPE(zsum(:,:,:), ish )
-              CALL mpp_sum(zwork, ish(1))
-              zsum(:,:,:)= RESHAPE(zwork,ish2)
-              DO jsec=1,nb_sec ; secs(jsec)%transport_h(:,:) = zsum(jsec,:,:) ; ENDDO
-           ENDIF
-
-           !Write the transport
-           DO jsec=1,nb_sec
-
-              IF( lwp .and.       ln_NOOS )CALL dia_dct_wri_NOOS_h(kt/nn_dctwri_h,jsec,secs(jsec))   ! use NOOS specific formatting
-            
-              !nullify transports values after writing
-              transports_3d_h(:,jsec,:,:)=0.0
-              transports_2d_h(:,jsec,:)=0.0
-              secs(jsec)%transport_h(:,:)=0.  
-              IF ( ln_NOOS ) CALL transport_h(secs(jsec),lldebug,jsec)  ! reinitialise for next 25 hour instantaneous average (overlapping values)
-
-           ENDDO
-
-        ENDIF 
-
-     ENDIF    
 
      IF( lk_mpp )THEN
@@ -419 +562 @@
      LOGICAL :: llbon                                       ,&!local logical
                 lldebug                                       !debug the section
+     LOGICAL       ::   verbose     
+     verbose = .false.
      !!-------------------------------------------------------------------------------------
      CALL wrk_alloc( nb_point_max, directemp )
@@ -424 +569 @@
      !open input file
      !---------------
-     CALL ctl_opn( numdct_in, 'section_ijglobal.diadct', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
+     !write(numout,*) 'dct low-level pre open: little endian '
+     !OPEN(UNIT=107,FILE='section_ijglobal.diadct', FORM='UNFORMATTED', ACCESS='SEQUENTIAL', STATUS='OLD',convert='LITTLE_ENDIAN')
+     
+     IF ( verbose ) write(numout,*) 'dct low-level pre open: big endian :',nproc,narea
+     OPEN(UNIT=107,FILE='section_ijglobal.diadct', FORM='UNFORMATTED', ACCESS='SEQUENTIAL', STATUS='OLD',convert='BIG_ENDIAN')
+     
+     !write(numout,*) 'dct low-level pre open: SWAP '
+     !OPEN(UNIT=107,FILE='section_ijglobal.diadct', FORM='UNFORMATTED', ACCESS='SEQUENTIAL', STATUS='OLD',convert='SWAP')
+     
+     !write(numout,*) 'dct low-level pre open: NATIVE '
+     !OPEN(UNIT=107,FILE='section_ijglobal.diadct', FORM='UNFORMATTED', ACCESS='SEQUENTIAL', STATUS='OLD',convert='NATIVE')
+     
+     READ(107) isec
+     CLOSE(107)
+     
+     CALL ctl_opn( numdct_in, 'section_ijglobal.diadct', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, .TRUE. )
  
      !---------------
@@ -433 +593 @@
 
         IF ( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) ) &
-           & WRITE(numout,*)'debuging for section number: ',jsec 
+           & WRITE(numout,*)'debugging for section number: ',jsec 
 
         !initialization
         !---------------
         secs(jsec)%name=''
-        secs(jsec)%llstrpond    = .FALSE.  ; secs(jsec)%ll_ice_section = .FALSE.
-        secs(jsec)%ll_date_line = .FALSE.  ; secs(jsec)%nb_class       = 0
-        secs(jsec)%zsigi        = 99._wp   ; secs(jsec)%zsigp          = 99._wp
-        secs(jsec)%zsal         = 99._wp   ; secs(jsec)%ztem           = 99._wp
-        secs(jsec)%zlay         = 99._wp         
-        secs(jsec)%transport    =  0._wp   ; secs(jsec)%nb_point       = 0
-        secs(jsec)%transport_h  =  0._wp   ; secs(jsec)%nb_point       = 0
+        secs(jsec)%llstrpond      = .FALSE.
+        secs(jsec)%ll_ice_section = .FALSE.
+        secs(jsec)%ll_date_line   = .FALSE.
+        secs(jsec)%nb_class       = 0
+        secs(jsec)%zsigi          = 99._wp
+        secs(jsec)%zsigp          = 99._wp
+        secs(jsec)%zsal           = 99._wp
+        secs(jsec)%ztem           = 99._wp
+        secs(jsec)%zlay           = 99._wp
+        secs(jsec)%transport      =  0._wp
+        secs(jsec)%transport_h    =  0._wp
+        secs(jsec)%nb_point       = 0
 
         !read section's number / name / computing choices / classes / slopeSection / points number
         !-----------------------------------------------------------------------------------------
-        READ(numdct_in,iostat=iost)isec
-        IF (iost .NE. 0 )EXIT !end of file
+        
+        READ(numdct_in,iostat=iost) isec
+        IF (iost .NE. 0 ) then
+          write(numout,*) 'reached end of section_ijglobal.diadct. iost = ',iost, &
+                          ', number of sections read = ', jsec-1
+          EXIT !end of file 
+        ENDIF
+        
+        
         WRITE(cltmp,'(a,i4.4,a,i4.4)')'diadct: read sections : Problem of section number: isec= ',isec,' and jsec= ',jsec
-        IF( jsec .NE. isec )CALL ctl_stop( cltmp )
-
-        IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )WRITE(numout,*)"isec ",isec 
+        
+        
+        IF( jsec .NE. isec )  CALL ctl_stop( cltmp )
 
         READ(numdct_in)secs(jsec)%name
@@ -469 +641 @@
         READ(numdct_in)iptglo
 
-        IF ( ln_NOOS ) THEN
+        IF ( ln_NOOS .AND. verbose ) THEN
            WRITE(numout,*) 'Section name = ',TRIM(secs(jsec)%name)
            WRITE(numout,*) 'coordSec = ',secs(jsec)%coordSec
@@ -483 +655 @@
 
             WRITE(numout,*)       "   Section name :                       ",TRIM(secs(jsec)%name)
-            WRITE(numout,*)       "      Compute heat and salt transport ? ",secs(jsec)%llstrpond
+            WRITE(numout,*)       "      Compute temperature and salinity transports ? ",secs(jsec)%llstrpond
             WRITE(numout,*)       "      Compute ice transport ?           ",secs(jsec)%ll_ice_section
             WRITE(numout,*)       "      Section crosses date-line ?       ",secs(jsec)%ll_date_line
@@ -500 +672 @@
            !read points'coordinates and directions 
            !--------------------------------------
-           IF ( ln_NOOS ) THEN
+           IF ( ln_NOOS .AND. verbose ) THEN
               WRITE(numout,*) 'Coords and direction read'
            ENDIF
@@ -558 +730 @@
            ENDIF
 
-              IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN
+           IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN
               DO jpt = 1,iptloc
                  iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
                  ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1
               ENDDO
-              ENDIF
+           ENDIF
 
            !remove redundant points between processors
@@ -602 +774 @@
  
      nb_sec = jsec-1   !number of section read in the file
+
+     IF( lwp .AND. verbose )  WRITE(numout,*)'diadct: read sections: Finished readsec.'
 
      CALL wrk_dealloc( nb_point_max, directemp )
@@ -703 +877 @@
      !!              loop on the level jk !!
      !! 
-     !!  Output  ::  Arrays containing the volume,density,heat,salt transports for each i
-     !!              point in a section, summed over each nn_dct. 
+     !! ** Output: Arrays containing the volume,density,salinity,temperature etc
+     !!            transports for each point in a section, summed over each nn_dct.
      !!
      !!-------------------------------------------------------------------------------------------
@@ -713 +887 @@
     
      !! * Local variables
-     INTEGER             :: jk, jseg, jclass,jl,                 &!loop on level/segment/classes/ice categories
-                            isgnu, isgnv                          ! 
-     REAL(wp)            :: zumid, zvmid,                        &!U/V velocity on a cell segment 
-                            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
+     INTEGER             :: jk,jseg,jclass,    &!loop on level/segment/classes 
+                            isgnu  , isgnv      !
+     REAL(wp):: zumid        , zvmid        ,  &!U/V velocity on a cell segment
+                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/ssh/potential density /depth at u/v point
 
      TYPE(POINT_SECTION) :: k
      !!--------------------------------------------------------
 
-     !!NIALL IF( ld_debug )WRITE(numout,*)'      Compute transport'
+     IF( ld_debug )WRITE(numout,*)'      Compute transport (jsec,sec%nb_point,sec%slopeSection) : ', jsec,sec%nb_point,sec%slopeSection
 
      !---------------------------!
@@ -730 +904 @@
      IF(sec%nb_point .NE. 0)THEN   
 
+        !----------------------------------------------------------------------------------------------------
+        !----------------------------------------------------------------------------------------------------
+        !----------------------------------------------------------------------------------------------------
+        !
+        !
+        ! ! ! ! JT 1/09/2018 - changing convention. Always direction + is toward left hand of section
+        !
+        !    Making sign of the velocities used to calculate the volume transport a function of direction, not slopesection
+        !    (isgnu, isgnv)
+        !    
+        !    They vary for each segment of the section. 
+        !
+        !----------------------------------------------------------------------------------------------------
+        !----------------------------------------------------------------------------------------------------
         !----------------------------------------------------------------------------------------------------
         !Compute sign for velocities:
@@ -751 +939 @@
         !                                                      
         !----------------------------------------------------------------------------------------------------
-        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
@@ -763 +946 @@
         !--------------------------------------!
         DO jseg=1,MAX(sec%nb_point-1,0)
-              
+           
+           
+           !Compute sign for velocities:
+           
+           !isgnu =  1
+           !isgnv =  1
+           !
+           !changing sign of u and v is dependent on the direction of the section. 
+           !isgnu =  1
+           !isgnv =  1
+           !SELECT CASE( sec%direction(jseg) )
+           !CASE(0)  ;   isgnv = -1
+           !CASE(3)  ;   isgnu = -1
+           !END SELECT
+           
+           
+           SELECT CASE( sec%direction(jseg) )
+           CASE(0)  
+              isgnu =  1
+              isgnv = -1
+           CASE(1)
+              isgnu =  1
+              isgnv =  1
+           CASE(2)  
+              isgnu =  1
+              isgnv =  1
+           CASE(3)  
+              isgnu = -1
+              isgnv =  1
+           END SELECT
+           
            !-------------------------------------------------------------------------------------------
            ! Select the appropriate coordinate for computing the velocity of the segment
+           ! Corrected by JT 01/09/2018 (#)
            !
            !                      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)
-           !                      --------                                              |
-           !                                                                            |
+           !      F(i,j)---------#V(i,j)-------F(i+1,j)                                 |
+           !                     -------->                                              |
+           !                                                                        |   |
+           !                                                                        |   |
+           !                      Case (3)                                          | U(i,j)
+           !                      --------                                          |   |
+           !                                                                        V   |
            !  listPoint(jseg+1) F(i,j+1)                                                |
            !                        |                                                   |
            !                        |                                                   |
            !                        |                                 listPoint(jseg+1) F(i,j-1)
-           !                        |                                            
-           !                        |                                            
-           !                     U(i,j+1)                                            
-           !                        |                                       Case(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)
+           !                        |                 F(i-1,j)----------#V(i-1,j) ------#f(i,j)                           
+           ! listPoint(jseg)     F(i,j)                                 <-------
            ! 
            !-------------------------------------------------------------------------------------------
@@ -800 +1014 @@
            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',0 ) 
-                 zsn   = interp(k%I,k%J,jk,'V',1 ) 
-                 zrhop = interp(k%I,k%J,jk,'V',2 ) 
-                 zrhoi = interp(k%I,k%J,jk,'V',3 ) 
-                 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',0 ) 
-                 zsn   = interp(k%I,k%J,jk,'U',1 ) 
-                 zrhop = interp(k%I,k%J,jk,'U',2 ) 
-                 zrhoi = interp(k%I,k%J,jk,'U',3 ) 
-                 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) 
-  
-              !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) 
+           !---------------------------|
+           !     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',0 )
+                 zsn   = interp(k%I,k%J,jk,'V',1 )
+                 zrhop = interp(k%I,k%J,jk,'V',2)
+                 zrhoi = interp(k%I,k%J,jk,'V',3)
+                 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',0)
+                 zsn   = interp(k%I,k%J,jk,'U',1)
+                 zrhop = interp(k%I,k%J,jk,'U',2)
+                 zrhoi = interp(k%I,k%J,jk,'U',3)
+                 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)
+ 
+              !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 
+              !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
               !COMPUTE TRANSPORT  
- 
-              transports_3d(1,jsec,jseg,jk) = transports_3d(1,jsec,jseg,jk) + zTnorm 
-  
-              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
+
+              transports_3d(1,jsec,jseg,jk) = transports_3d(1,jsec,jseg,jk) + zTnorm
+ 
+              IF ( sec%llstrpond ) THEN
+                 transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk)  + zTnorm * zrhoi
+                 transports_3d(3,jsec,jseg,jk) = transports_3d(3,jsec,jseg,jk)  + zTnorm * zrhop
                  transports_3d(4,jsec,jseg,jk) = transports_3d(4,jsec,jseg,jk)  + zTnorm * 4.e+3_wp * (ztn+273.15) * 1026._wp
                  transports_3d(5,jsec,jseg,jk) = transports_3d(5,jsec,jseg,jk)  + zTnorm * 0.001 * zsn * 1026._wp
@@ -880 +1094 @@
    
               zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J)
-
-#if defined key_lim2   
-              transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)*   & 
-                                   (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  & 
-                                  *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) +  & 
-                                    hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))
-              transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)*   & 
-                                    (1.0 -  frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))
-#endif
-#if defined key_lim3
-              DO jl=1,jpl
-                 transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)*     &
-                                   a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) * &
-                                  ( ht_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) +  &
-                                    ht_s(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) )
-                                   
-                 transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)*   &
-                                   a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl)
-              ENDDO
-#endif
+   
+              transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)*   &
+                                   (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &
+                                  *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) +  &
+                                    hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) &
+                                   +zice_vol_pos
+              transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)*   &
+                                    (1.0 -  frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &
+                                   +zice_surf_pos
    
            ENDIF !end of ice case
@@ -957 +1160 @@
 
         !----------------------------------------------------------------------------------------------------
+        !----------------------------------------------------------------------------------------------------
+        !----------------------------------------------------------------------------------------------------
+        !
+        !
+        ! ! ! ! JT 1/09/2018 - changing convention. Always direction + is toward left hand of section
+        !
+        !    Making sign of the velocities used to calculate the volume transport a function of direction, not slopesection
+        !    (isgnu, isgnv)
+        !    
+        !    They vary for each segment of the section. 
+        !
+        !----------------------------------------------------------------------------------------------------
+        !----------------------------------------------------------------------------------------------------
+        !----------------------------------------------------------------------------------------------------
         !Compute sign for velocities:
         !
@@ -977 +1194 @@
         !                                                      
         !----------------------------------------------------------------------------------------------------
-        isgnu = 1
-        IF( sec%slopeSection .GT. 0 ) THEN  ; isgnv = -1 
-        ELSE                                ; isgnv =  1
-        ENDIF
 
         IF( ld_debug )write(numout,*)"isgnu isgnv ",isgnu,isgnv
@@ -988 +1201 @@
         !--------------------------------------!
         DO jseg=1,MAX(sec%nb_point-1,0)
-              
+           
+           
+           !Compute sign for velocities:
+           
+           !isgnu =  1
+           !isgnv =  1
+           !
+           ! changing sign of u and v is dependent on the direction of the section. 
+           !isgnu =  1
+           !isgnv =  1
+           !SELECT CASE( sec%direction(jseg) )
+           !CASE(0)  ;   isgnv = -1
+           !CASE(3)  ;   isgnu = -1
+           !END SELECT
+           
+           
+           SELECT CASE( sec%direction(jseg) )
+           CASE(0)  
+              isgnu =  1
+              isgnv = -1
+           CASE(1)
+              isgnu =  1
+              isgnv =  1
+           CASE(2)  
+              isgnu =  1
+              isgnv =  1
+           CASE(3)  
+              isgnu = -1
+              isgnv =  1
+           END SELECT
+           
            !-------------------------------------------------------------------------------------------
            ! Select the appropriate coordinate for computing the velocity of the segment
+           ! Corrected by JT 01/09/2018 (#)
            !
            !                      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)
-           !                      --------                                              |
-           !                                                                            |
+           !      F(i,j)---------#V(i,j)-------F(i+1,j)                                 |
+           !                     -------->                                              |
+           !                                                                        |   |
+           !                                                                        |   |
+           !                      Case (3)                                          | U(i,j)
+           !                      --------                                          |   |
+           !                                                                        V   |
            !  listPoint(jseg+1) F(i,j+1)                                                |
            !                        |                                                   |
            !                        |                                                   |
            !                        |                                 listPoint(jseg+1) F(i,j-1)
-           !                        |                                            
-           !                        |                                            
-           !                     U(i,j+1)                                            
-           !                        |                                       Case(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)
+           !                        |                 F(i-1,j)----------#V(i-1,j) ------#f(i,j)                           
+           ! listPoint(jseg)     F(i,j)                                 <-------
            ! 
            !-------------------------------------------------------------------------------------------
@@ -1031 +1275 @@
            DO jk=1,mbathy(k%I,k%J)
 
-              ! compute temparature, salinity, insitu & potential density, ssh and depth at U/V point
+              ! 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',0 )
+                 ztn   = interp(k%I,k%J,jk,'V',0)
                  zsn   = interp(k%I,k%J,jk,'V',1)
                  zrhop = interp(k%I,k%J,jk,'V',2)
@@ -1151 +1395 @@
      REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point 
      !!------------------------------------------------------------- 
- 
-     !! Sum the relevant segments to obtain values for each class 
-     IF(sec%nb_point .NE. 0)THEN    
- 
-        !--------------------------------------! 
-        ! 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',0 ) 
-                 zsn   = interp(k%I,k%J,jk,'V',1 ) 
-                 zrhop = interp(k%I,k%J,jk,'V',2) 
-                 zrhoi = interp(k%I,k%J,jk,'V',3) 
-
-              CASE(2,3) 
-                 ztn   = interp(k%I,k%J,jk,'U',0 ) 
-                 zsn   = interp(k%I,k%J,jk,'U',1 ) 
-                 zrhop = interp(k%I,k%J,jk,'U',2 ) 
-                 zrhoi = interp(k%I,k%J,jk,'U',3 ) 
-                 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) 
-  
-              !------------------------------- 
-              !  LOOP ON THE DENSITY CLASSES | 
-              !------------------------------- 
-              !The computation is made for each density/temperature/salinity/depth class 
-              DO jclass=1,MAX(1,sec%nb_class-1) 
- 
-                 !----------------------------------------------! 
-                 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL!  
-                 !----------------------------------------------! 
-
-                 IF ( (                                                    & 
-                    ((( zrhop .GE. (sec%zsigp(jclass)+1000.  )) .AND.      & 
-                    (   zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR.     & 
-                    ( sec%zsigp(jclass) .EQ. 99.)) .AND.                   & 
- 
-                    ((( zrhoi .GE. (sec%zsigi(jclass) + 1000.  )) .AND.    & 
-                    (   zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR.     & 
-                    ( sec%zsigi(jclass) .EQ. 99.)) .AND.                   & 
- 
-                    ((( zsn .GT. sec%zsal(jclass)) .AND.                   & 
-                    (   zsn .LE. sec%zsal(jclass+1))) .OR.                 & 
-                    ( sec%zsal(jclass) .EQ. 99.)) .AND.                    & 
- 
-                    ((( ztn .GE. sec%ztem(jclass)) .AND.                   & 
-                    (   ztn .LE. sec%ztem(jclass+1))) .OR.                 & 
-                    ( sec%ztem(jclass) .EQ.99.)) .AND.                     & 
- 
-                    ((( zfsdep .GE. sec%zlay(jclass)) .AND.                & 
-                    (   zfsdep .LE. sec%zlay(jclass+1))) .OR.              & 
-                    ( sec%zlay(jclass) .EQ. 99. ))                         & 
-                                                                   ))   THEN 
- 
-                    !SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS 
-                    !---------------------------------------------------------------------------- 
-                    IF (transports_3d(1,jsec,jseg,jk) .GE. 0.0) THEN  
-                       sec%transport(1,jclass) = sec%transport(1,jclass)+transports_3d(1,jsec,jseg,jk)*1.E-6 
-                    ELSE 
-                       sec%transport(2,jclass) = sec%transport(2,jclass)+transports_3d(1,jsec,jseg,jk)*1.E-6 
-                    ENDIF 
-                    IF( sec%llstrpond )THEN 
- 
-                       IF ( transports_3d(2,jsec,jseg,jk) .GE. 0.0 ) THEN 
-                          sec%transport(3,jclass) = sec%transport(3,jclass)+transports_3d(2,jsec,jseg,jk) 
-                       ELSE 
-                          sec%transport(4,jclass) = sec%transport(4,jclass)+transports_3d(2,jsec,jseg,jk) 
-                       ENDIF 
- 
-                       IF ( transports_3d(3,jsec,jseg,jk) .GE. 0.0 ) THEN 
-                          sec%transport(5,jclass) = sec%transport(5,jclass)+transports_3d(3,jsec,jseg,jk) 
-                       ELSE 
-                          sec%transport(6,jclass) = sec%transport(6,jclass)+transports_3d(3,jsec,jseg,jk) 
-                       ENDIF 
-
-                       IF ( transports_3d(4,jsec,jseg,jk) .GE. 0.0 ) THEN
-                          sec%transport(7,jclass) = sec%transport(7,jclass)+transports_3d(4,jsec,jseg,jk)
-                       ELSE
-                          sec%transport(8,jclass) = sec%transport(8,jclass)+transports_3d(4,jsec,jseg,jk)
-                       ENDIF
-
-                       IF ( transports_3d(5,jsec,jseg,jk) .GE. 0.0 ) THEN
-                          sec%transport( 9,jclass) = sec%transport( 9,jclass)+transports_3d(5,jsec,jseg,jk)
-                       ELSE
-                          sec%transport(10,jclass) = sec%transport(10,jclass)+transports_3d(5,jsec,jseg,jk)
-                       ENDIF
- 
-                    ELSE 
-                       sec%transport( 3,jclass) = 0._wp 
-                       sec%transport( 4,jclass) = 0._wp 
-                       sec%transport( 5,jclass) = 0._wp 
-                       sec%transport( 6,jclass) = 0._wp 
-                       sec%transport( 7,jclass) = 0._wp
-                       sec%transport( 8,jclass) = 0._wp
-                       sec%transport( 9,jclass) = 0._wp
-                       sec%transport(10,jclass) = 0._wp
-                    ENDIF 
- 
-                 ENDIF ! end of test if point is in class 
-    
-              ENDDO ! end of loop on the classes 
- 
-           ENDDO ! loop over jk 
- 
-#if defined key_lim2 || defined key_lim3 
- 
-           !ICE CASE     
-           IF( sec%ll_ice_section )THEN 
- 
-              IF ( transports_2d(1,jsec,jseg) .GE. 0.0 ) THEN 
-                 sec%transport(11,1) = sec%transport(11,1)+transports_2d(1,jsec,jseg)*1.E-6 
-              ELSE 
-                 sec%transport(12,1) = sec%transport(12,1)+transports_2d(1,jsec,jseg)*1.E-6 
-              ENDIF 
- 
-              IF ( transports_2d(3,jsec,jseg) .GE. 0.0 ) THEN 
-                 sec%transport(13,1) = sec%transport(13,1)+transports_2d(2,jsec,jseg)*1.E-6 
-              ELSE 
-                 sec%transport(14,1) = sec%transport(14,1)+transports_2d(2,jsec,jseg)*1.E-6 
-              ENDIF 
- 
-           ENDIF !end of ice case 
-#endif 
-  
-        ENDDO !end of loop on the segment 
- 
-     ELSE  !if sec%nb_point =0 
-        sec%transport(1:2,:)=0. 
-        IF (sec%llstrpond) sec%transport(3:10,:)=0. 
-        IF (sec%ll_ice_section) sec%transport(11:14,:)=0. 
-     ENDIF !end of sec%nb_point =0 case 
- 
-  END SUBROUTINE dia_dct_sum 
-
-  SUBROUTINE dia_dct_sum_h(sec,jsec)
-     !!-------------------------------------------------------------
-     !! Exactly as dia_dct_sum but for hourly files containing data summed at each time step
-     !!
-     !! Purpose: Average the transport over nn_dctwri time steps 
-     !! and sum over the density/salinity/temperature/depth classes
-     !!
-     !! Method: 
-     !!           Sum over relevant grid cells to obtain values
-     !!           for each
-     !!              There are several loops:                 
-     !!              loop on the segment between 2 nodes
-     !!              loop on the level jk
-     !!              loop on the density/temperature/salinity/level classes
-     !!              test on the density/temperature/salinity/level
-     !!
-     !!  ** Method  :Transport through a given section is equal to the sum of transports
-     !!              computed on each proc.
-     !!              On each proc,transport is equal to the sum of transport computed through
-     !!              segments linking each point of sec%listPoint  with the next one.   
-     !!
-     !!-------------------------------------------------------------
-     !! * arguments
-     TYPE(SECTION),INTENT(INOUT) :: sec
-     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section
-
-     TYPE(POINT_SECTION) :: k
-     INTEGER  :: jk,jseg,jclass                        !loop on level/segment/classes 
-     REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point
-     !!-------------------------------------------------------------
+
 
      !! Sum the relevant segments to obtain values for each class
@@ -1416 +1450 @@
               SELECT CASE( sec%direction(jseg) )
               CASE(0,1)
-                 ztn   = interp(k%I,k%J,jk,'V',0 )
-                 zsn   = interp(k%I,k%J,jk,'V',1 )
-                 zrhop = interp(k%I,k%J,jk,'V',2 )
-                 zrhoi = interp(k%I,k%J,jk,'V',3 )
+                 ztn   = interp(k%I,k%J,jk,'V',0)
+                 zsn   = interp(k%I,k%J,jk,'V',1)
+                 zrhop = interp(k%I,k%J,jk,'V',2)
+                 zrhoi = interp(k%I,k%J,jk,'V',3)
                  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',0 )
-                 zsn   = interp(k%I,k%J,jk,'U',1 )
-                 zrhop = interp(k%I,k%J,jk,'U',2 )
-                 zrhoi = interp(k%I,k%J,jk,'U',3 )
+                 ztn   = interp(k%I,k%J,jk,'U',0)
+                 zsn   = interp(k%I,k%J,jk,'U',1)
+                 zrhop = interp(k%I,k%J,jk,'U',2)
+                 zrhoi = interp(k%I,k%J,jk,'U',3)
+                 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) 
+ 
+              !-------------------------------
+              !  LOOP ON THE DENSITY CLASSES |
+              !-------------------------------
+              !The computation is made for each density/temperature/salinity/depth class
+              DO jclass=1,MAX(1,sec%nb_class-1)
+ 
+                 !----------------------------------------------!
+                 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL! 
+                 !----------------------------------------------!
+ 
+                 IF ( (                                                    &
+                    ((( zrhop .GE. (sec%zsigp(jclass)+1000.  )) .AND.      &
+                    (   zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR.     &
+                    ( sec%zsigp(jclass) .EQ. 99.)) .AND.                   &
+
+                    ((( zrhoi .GE. (sec%zsigi(jclass) + 1000.  )) .AND.    &
+                    (   zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR.     &
+                    ( sec%zsigi(jclass) .EQ. 99.)) .AND.                   &
+
+                    ((( zsn .GT. sec%zsal(jclass)) .AND.                   &
+                    (   zsn .LE. sec%zsal(jclass+1))) .OR.                 &
+                    ( sec%zsal(jclass) .EQ. 99.)) .AND.                    &
+
+                    ((( ztn .GE. sec%ztem(jclass)) .AND.                   &
+                    (   ztn .LE. sec%ztem(jclass+1))) .OR.                 &
+                    ( sec%ztem(jclass) .EQ.99.)) .AND.                     &
+
+                    ((( zfsdep .GE. sec%zlay(jclass)) .AND.                &
+                    (   zfsdep .LE. sec%zlay(jclass+1))) .OR.              &
+                    ( sec%zlay(jclass) .EQ. 99. ))                         &
+                                                                   ))   THEN
+
+                    !SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS
+                    !----------------------------------------------------------------------------
+                    IF (transports_3d(1,jsec,jseg,jk) .GE. 0.0) THEN 
+                       sec%transport(1,jclass) = sec%transport(1,jclass)+transports_3d(1,jsec,jseg,jk)
+                    ELSE
+                       sec%transport(2,jclass) = sec%transport(2,jclass)+transports_3d(1,jsec,jseg,jk)
+                    ENDIF
+                    IF( sec%llstrpond )THEN
+
+                       IF( transports_3d(1,jsec,jseg,jk) .NE. 0._wp ) THEN
+
+                          IF (transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) .GE. 0.0 ) THEN
+                             sec%transport(3,jclass) = sec%transport(3,jclass)+transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk)
+                          ELSE
+                             sec%transport(4,jclass) = sec%transport(4,jclass)+transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk)
+                          ENDIF
+
+                          IF ( transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) .GE. 0.0 ) THEN
+                             sec%transport(5,jclass) = sec%transport(5,jclass)+transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk)
+                          ELSE
+                             sec%transport(6,jclass) = sec%transport(6,jclass)+transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk)
+                          ENDIF
+
+                       ENDIF
+
+                       IF ( transports_3d(4,jsec,jseg,jk) .GE. 0.0 ) THEN
+                          sec%transport(7,jclass) = sec%transport(7,jclass)+transports_3d(4,jsec,jseg,jk)
+                       ELSE
+                          sec%transport(8,jclass) = sec%transport(8,jclass)+transports_3d(4,jsec,jseg,jk)
+                       ENDIF
+
+                       IF ( transports_3d(5,jsec,jseg,jk) .GE. 0.0 ) THEN
+                          sec%transport( 9,jclass) = sec%transport( 9,jclass)+transports_3d(5,jsec,jseg,jk)
+                       ELSE
+                          sec%transport(10,jclass) = sec%transport(10,jclass)+transports_3d(5,jsec,jseg,jk)
+                       ENDIF
+ 
+                    ELSE 
+                       sec%transport( 3,jclass) = 0._wp 
+                       sec%transport( 4,jclass) = 0._wp 
+                       sec%transport( 5,jclass) = 0._wp 
+                       sec%transport( 6,jclass) = 0._wp 
+                       sec%transport( 7,jclass) = 0._wp
+                       sec%transport( 8,jclass) = 0._wp
+                       sec%transport( 9,jclass) = 0._wp
+                       sec%transport(10,jclass) = 0._wp
+                    ENDIF 
+ 
+                 ENDIF ! end of test if point is in class 
+    
+              ENDDO ! end of loop on the classes 
+ 
+           ENDDO ! loop over jk 
+ 
+#if defined key_lim2 || defined key_lim3 
+ 
+           !ICE CASE     
+           IF( sec%ll_ice_section )THEN 
+ 
+              IF ( transports_2d(1,jsec,jseg) .GE. 0.0 ) THEN 
+                 sec%transport(11,1) = sec%transport(11,1)+transports_2d(1,jsec,jseg)
+              ELSE 
+                 sec%transport(12,1) = sec%transport(12,1)+transports_2d(1,jsec,jseg)
+              ENDIF 
+ 
+              IF ( transports_2d(3,jsec,jseg) .GE. 0.0 ) THEN 
+                 sec%transport(13,1) = sec%transport(13,1)+transports_2d(2,jsec,jseg)
+              ELSE 
+                 sec%transport(14,1) = sec%transport(14,1)+transports_2d(2,jsec,jseg)
+              ENDIF 
+ 
+           ENDIF !end of ice case 
+#endif 
+ 
+        ENDDO !end of loop on the segment
+
+     ELSE  !if sec%nb_point =0
+        sec%transport(1:2,:)=0.
+        IF (sec%llstrpond) sec%transport(3:10,:)=0.
+        IF (sec%ll_ice_section) sec%transport( 11:14,:)=0.
+     ENDIF !end of sec%nb_point =0 case
+
+  END SUBROUTINE dia_dct_sum
+ 
+  SUBROUTINE dia_dct_sum_h(sec,jsec)
+     !!-------------------------------------------------------------
+     !! Exactly as dia_dct_sum but for hourly files containing data summed at each time step
+     !!
+     !! Purpose: Average the transport over nn_dctwri time steps 
+     !! and sum over the density/salinity/temperature/depth classes
+     !!
+     !! Method: 
+     !!           Sum over relevant grid cells to obtain values
+     !!           for each
+     !!              There are several loops:                 
+     !!              loop on the segment between 2 nodes
+     !!              loop on the level jk
+     !!              loop on the density/temperature/salinity/level classes
+     !!              test on the density/temperature/salinity/level
+     !!
+     !!  ** Method  :Transport through a given section is equal to the sum of transports
+     !!              computed on each proc.
+     !!              On each proc,transport is equal to the sum of transport computed through
+     !!              segments linking each point of sec%listPoint  with the next one.   
+     !!
+     !!-------------------------------------------------------------
+     !! * arguments
+     TYPE(SECTION),INTENT(INOUT) :: sec
+     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section
+
+     TYPE(POINT_SECTION) :: k
+     INTEGER  :: jk,jseg,jclass                        !loop on level/segment/classes 
+     REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point
+     !!-------------------------------------------------------------
+
+     !! Sum the relevant segments to obtain values for each class
+     IF(sec%nb_point .NE. 0)THEN   
+
+        !--------------------------------------!
+        ! 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',0)
+                 zsn   = interp(k%I,k%J,jk,'V',1)
+                 zrhop = interp(k%I,k%J,jk,'V',2)
+                 zrhoi = interp(k%I,k%J,jk,'V',3)
+                 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',0)
+                 zsn   = interp(k%I,k%J,jk,'U',1)
+                 zrhop = interp(k%I,k%J,jk,'U',2)
+                 zrhoi = interp(k%I,k%J,jk,'U',3)
                  zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1) 
               END SELECT
@@ -1472 +1721 @@
                     IF( sec%llstrpond )THEN
 
-                       IF ( transports_3d_h(2,jsec,jseg,jk) .GE. 0.0 ) THEN 
-                          sec%transport_h(3,jclass) = sec%transport_h(3,jclass)+transports_3d_h(2,jsec,jseg,jk) 
-                       ELSE 
-                          sec%transport_h(4,jclass) = sec%transport_h(4,jclass)+transports_3d_h(2,jsec,jseg,jk) 
-                       ENDIF 
- 
-                       IF ( transports_3d_h(3,jsec,jseg,jk) .GE. 0.0 ) THEN 
-                          sec%transport_h(5,jclass) = sec%transport_h(5,jclass)+transports_3d_h(3,jsec,jseg,jk) 
-                       ELSE 
-                          sec%transport_h(6,jclass) = sec%transport_h(6,jclass)+transports_3d_h(3,jsec,jseg,jk) 
-                       ENDIF 
+                       IF( transports_3d_h(1,jsec,jseg,jk) .NE. 0._wp ) THEN
+
+                          IF (transports_3d_h(2,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) .GE. 0.0 ) THEN
+                             sec%transport_h(3,jclass) = sec%transport_h(3,jclass)+transports_3d_h(2,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk)
+                          ELSE
+                             sec%transport_h(4,jclass) = sec%transport_h(4,jclass)+transports_3d_h(2,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk)
+                          ENDIF
+
+                          IF ( transports_3d_h(3,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) .GE. 0.0 ) THEN
+                             sec%transport_h(5,jclass) = sec%transport_h(5,jclass)+transports_3d_h(3,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk)
+                          ELSE
+                             sec%transport_h(6,jclass) = sec%transport_h(6,jclass)+transports_3d_h(3,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk)
+                          ENDIF
+
+                       ENDIF
 
                        IF ( transports_3d_h(4,jsec,jseg,jk) .GE. 0.0 ) THEN
@@ -1572 +1825 @@
      !
      REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace 
+     CHARACTER(len=3)      :: noos_sect_name             ! Classname 
+     CHARACTER(len=25)      :: noos_var_sect_name
+     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   noos_iom_dummy
+     INTEGER               :: IERR
+     
+     REAL(wp), DIMENSION(3) :: tmp_iom_output
+     REAL(wp)               :: max_iom_val
+     LOGICAL       ::   verbose     
+     verbose = .false.
+     
      !!------------------------------------------------------------- 
+     
+     
+     
+     IF( lwp .AND. verbose ) THEN
+        WRITE(numout,*) " "
+        WRITE(numout,*) "dia_dct_wri_NOOS: write transports through sections at timestep: ", kt
+        WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
+     ENDIF   
+        
      CALL wrk_alloc(nb_type_class , zsumclasses )  
 
      zsumclasses(:)=0._wp
      zslope = sec%slopeSection       
-
-     WRITE(numdct_NOOS,'(I4,a1,I2,a1,I2,a12,i3,a17,i3,a10,a25)') nyear,'.',nmonth,'.',nday,'   Transect:',ksec-1,'   No. of layers:',sec%nb_class-1,'   Name: ',sec%name
-
+     
+     IF( lwp ) THEN
+         IF  ( ln_dct_ascii ) THEN
+             WRITE(numdct_NOOS,'(I4,a1,I2,a1,I2,a12,i3,a17,i3,a10,a25)') nyear,'.',nmonth,'.',nday,'   Transect:',ksec-1,'   No. of layers:',sec%nb_class-1,'   Name: ',sec%name
+         ELSE
+             WRITE(numdct_NOOS) nyear,nmonth,nday,ksec-1,sec%nb_class-1,sec%name
+         ENDIF  
+     ENDIF
+    
+     ! Sum all classes together, to give one values per type (pos tran, neg vol trans etc...). 
      DO jclass=1,MAX(1,sec%nb_class-1)
         zsumclasses(1:nb_type_class)=zsumclasses(1:nb_type_class)+sec%transport(1:nb_type_class,jclass)
@@ -1587 +1866 @@
      zbnd1   = 0._wp
      zbnd2   = 0._wp
-
-     IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN
-        WRITE(numdct_NOOS,'(9e12.4E2)') -(zsumclasses( 1)+zsumclasses( 2)), -zsumclasses( 2),-zsumclasses( 1),   &
-                                        -(zsumclasses( 7)+zsumclasses( 8)), -zsumclasses( 8),-zsumclasses( 7),   &
-                                        -(zsumclasses( 9)+zsumclasses(10)), -zsumclasses(10),-zsumclasses( 9)
-     ELSE
-        WRITE(numdct_NOOS,'(9e12.4E2)')   zsumclasses( 1)+zsumclasses( 2) ,  zsumclasses( 1), zsumclasses( 2),   &
-                                          zsumclasses( 7)+zsumclasses( 8) ,  zsumclasses( 7), zsumclasses( 8),   &
-                                          zsumclasses( 9)+zsumclasses(10) ,  zsumclasses( 9), zsumclasses(10)
-     ENDIF 
+     
+     
+     
+     write (noos_sect_name, "(I0.3)")  ksec
+     
+     IF ( ln_dct_iom_cont ) THEN
+         max_iom_val = 1.e10
+         ALLOCATE( noos_iom_dummy(jpi,jpj,3),  STAT= ierr )
+            IF( ierr /= 0 )   CALL ctl_stop( 'dia_dct_wri_NOOS: failed to allocate noos_iom_dummy array' )
+     ENDIF
+     
+!     JT   I think changing the sign on the output based on the zslope value is redunant.
+!
+!     IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN
+!         
+!         IF( lwp ) THEN
+!             WRITE(numdct_NOOS,'(9e12.4E2)') -(zsumclasses( 1)+zsumclasses( 2)), -zsumclasses( 2),-zsumclasses( 1),   &
+!                                            -(zsumclasses( 7)+zsumclasses( 8)), -zsumclasses( 8),-zsumclasses( 7),   &
+!                                            -(zsumclasses( 9)+zsumclasses(10)), -zsumclasses(10),-zsumclasses( 9)
+!             CALL FLUSH(numdct_NOOS) 
+!         endif
+
+!         
+!         IF ( ln_dct_iom_cont ) THEN
+!         
+!             noos_iom_dummy(:,:,:) = 0.
+!             
+!             tmp_iom_output(:) = 0.
+!             tmp_iom_output(1) = -(zsumclasses( 1)+zsumclasses( 2))
+!             tmp_iom_output(2) = -zsumclasses( 2)
+!             tmp_iom_output(3) = -zsumclasses( 1)
+!             
+!             ! Convert to Sv 
+!             tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
+!             tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
+!             tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
+!             
+!             ! limit maximum and minimum values in iom_put
+!             if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
+!             if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
+!             if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
+!             if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
+!             if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
+!             if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
+!             
+!             ! Set NaN's to Zero         
+!             if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
+!             if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
+!             if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
+!             
+!             noos_iom_dummy(:,:,1) = tmp_iom_output(1)
+!             noos_iom_dummy(:,:,2) = tmp_iom_output(2)
+!             noos_iom_dummy(:,:,3) = tmp_iom_output(3)
+!             
+!             noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_trans'
+!             if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
+!             CALL iom_put( noos_var_sect_name,  noos_iom_dummy )
+!             noos_iom_dummy(:,:,:) = 0.         
+!             tmp_iom_output(:) = 0.
+!             tmp_iom_output(1) = -(zsumclasses( 7)+zsumclasses( 8))
+!             tmp_iom_output(2) = -zsumclasses( 8)
+!             tmp_iom_output(3) = -zsumclasses( 7)
+!             
+!             ! Convert to TJ/s
+!             tmp_iom_output(1) = tmp_iom_output(1)*1.E-12
+!             tmp_iom_output(2) = tmp_iom_output(2)*1.E-12
+!             tmp_iom_output(3) = tmp_iom_output(3)*1.E-12
+!             
+!             ! limit maximum and minimum values in iom_put
+!             if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
+!             if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
+!             if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
+!             if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
+!             if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
+!             if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
+!             
+!             ! Set NaN's to Zero         
+!             if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
+!             if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
+!             if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
+!             
+!             noos_iom_dummy(:,:,1) = tmp_iom_output(1)
+!             noos_iom_dummy(:,:,2) = tmp_iom_output(2)
+!             noos_iom_dummy(:,:,3) = tmp_iom_output(3)
+!             
+!             !noos_iom_dummy(:,:,1) = -(zsumclasses( 7)+zsumclasses( 8))
+!             !noos_iom_dummy(:,:,2) = -zsumclasses( 8)
+!             !noos_iom_dummy(:,:,3) = -zsumclasses( 7)
+!             
+!             noos_var_sect_name =  "noos_" // trim(noos_sect_name) // '_heat'
+!             if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
+!             CALL iom_put( noos_var_sect_name,  noos_iom_dummy )
+!              
+!             noos_iom_dummy(:,:,:) = 0.         
+!             tmp_iom_output(:) = 0.         
+!             tmp_iom_output(1) = -(zsumclasses( 9)+zsumclasses( 10))
+!             tmp_iom_output(2) = -zsumclasses( 10)
+!             tmp_iom_output(3) = -zsumclasses( 9)
+!             
+!             ! Convert to  MT/s
+!             tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
+!             tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
+!             tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
+!             
+!             ! limit maximum and minimum values in iom_put
+!             if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
+!             if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
+!             if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
+!             if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
+!             if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
+!             if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
+!             
+!             ! Set NaN's to Zero         
+!             if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
+!             if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
+!             if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
+!             
+!             noos_iom_dummy(:,:,1) = tmp_iom_output(1)
+!             noos_iom_dummy(:,:,2) = tmp_iom_output(2)
+!             noos_iom_dummy(:,:,3) = tmp_iom_output(3)
+!             
+!             !noos_iom_dummy(:,:,1) = -(zsumclasses( 9)+zsumclasses( 10))
+!             !noos_iom_dummy(:,:,2) = -zsumclasses( 10)
+!             !noos_iom_dummy(:,:,3) = -zsumclasses( 9)
+!             
+!             noos_var_sect_name =  "noos_" // trim(noos_sect_name) // '_salt'
+!             if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
+!             CALL iom_put( noos_var_sect_name,  noos_iom_dummy )
+!             noos_iom_dummy(:,:,:) = 0.     
+!             tmp_iom_output(:) = 0.    
+!         ENDIF
+!     ELSE
+!        IF( lwp ) THEN
+!            WRITE(numdct_NOOS,'(9e12.4E2)')   zsumclasses( 1)+zsumclasses( 2) ,  zsumclasses( 1), zsumclasses( 2),   &
+!                                            zsumclasses( 7)+zsumclasses( 8) ,  zsumclasses( 7), zsumclasses( 8),   &
+!                                            zsumclasses( 9)+zsumclasses(10) ,  zsumclasses( 9), zsumclasses(10)
+!            CALL FLUSH(numdct_NOOS) 
+!        endif
+!         
+!         
+!        IF ( ln_dct_iom_cont ) THEN
+!            
+!             noos_iom_dummy(:,:,:) = 0.
+!             tmp_iom_output(:) = 0.
+!             
+!             tmp_iom_output(1) = (zsumclasses( 1)+zsumclasses( 2))
+!             tmp_iom_output(2) = zsumclasses( 1)
+!             tmp_iom_output(3) = zsumclasses( 2)
+!             
+!             ! Convert to Sv
+!             tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
+!             tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
+!             tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
+!             
+!             ! limit maximum and minimum values in iom_put
+!             if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
+!             if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
+!             if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
+!             if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
+!             if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
+!             if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
+!             
+!             ! Set NaN's to Zero         
+!             if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
+!             if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
+!             if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
+!             
+!             noos_iom_dummy(:,:,1) = tmp_iom_output(1)
+!             noos_iom_dummy(:,:,2) = tmp_iom_output(2)
+!             noos_iom_dummy(:,:,3) = tmp_iom_output(3)
+!             
+!             !noos_iom_dummy(:,:,1) = (zsumclasses( 1)+zsumclasses( 2))
+!             !noos_iom_dummy(:,:,2) = zsumclasses( 1)
+!             !noos_iom_dummy(:,:,3) = zsumclasses( 2)
+!             
+!             
+!             
+!             noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_trans'
+!             if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
+!             CALL iom_put( noos_var_sect_name,  noos_iom_dummy )
+!             noos_iom_dummy(:,:,:) = 0.
+!             tmp_iom_output(:) = 0.
+!             
+!             tmp_iom_output(1) = (zsumclasses( 7)+zsumclasses( 8))
+!             tmp_iom_output(2) = zsumclasses( 7)
+!             tmp_iom_output(3) = zsumclasses( 8)
+!             
+!             ! Convert to TJ/s
+!             tmp_iom_output(1) = tmp_iom_output(1)*1.E-12
+!             tmp_iom_output(2) = tmp_iom_output(2)*1.E-12
+!             tmp_iom_output(3) = tmp_iom_output(3)*1.E-12
+!             
+!             ! limit maximum and minimum values in iom_put
+!             if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
+!             if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
+!             if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
+!             if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
+!             if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
+!             if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
+!             
+!             ! Set NaN's to Zero         
+!             if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
+!             if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
+!             if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
+!             
+!             noos_iom_dummy(:,:,1) = tmp_iom_output(1)
+!             noos_iom_dummy(:,:,2) = tmp_iom_output(2)
+!             noos_iom_dummy(:,:,3) = tmp_iom_output(3)
+!             
+!             !noos_iom_dummy(:,:,1) = (zsumclasses( 7)+zsumclasses( 8))
+!             !noos_iom_dummy(:,:,2) = zsumclasses( 7)
+!             !noos_iom_dummy(:,:,3) = zsumclasses( 8)
+!             
+!             noos_var_sect_name =  "noos_" // trim(noos_sect_name) // '_heat'
+!             if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
+!             CALL iom_put(noos_var_sect_name,  noos_iom_dummy )  
+!             noos_iom_dummy(:,:,:) = 0.
+!             tmp_iom_output(:) = 0.
+!             
+!             tmp_iom_output(1) = (zsumclasses( 9)+zsumclasses( 10))
+!             tmp_iom_output(2) = zsumclasses( 9)
+!             tmp_iom_output(3) = zsumclasses( 10)
+!             
+!             ! Convert to  MT/s
+!             tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
+!             tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
+!             tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
+!             
+!             
+!             ! limit maximum and minimum values in iom_put
+!             if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
+!             if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
+!             if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
+!             if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
+!             if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
+!             if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
+!             
+!             ! Set NaN's to Zero         
+!             if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
+!             if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
+!             if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
+!             
+!             noos_iom_dummy(:,:,1) = tmp_iom_output(1)
+!             noos_iom_dummy(:,:,2) = tmp_iom_output(2)
+!             noos_iom_dummy(:,:,3) = tmp_iom_output(3)
+!             
+!             !noos_iom_dummy(:,:,1) = (zsumclasses( 9)+zsumclasses( 10))
+!             !noos_iom_dummy(:,:,2) = zsumclasses( 9)
+!             !noos_iom_dummy(:,:,3) = zsumclasses( 10)
+!             
+!             noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_salt'
+!             if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
+!             CALL iom_put(noos_var_sect_name,  noos_iom_dummy )
+!             noos_iom_dummy(:,:,:) = 0.         
+!             tmp_iom_output(:) = 0.
+!         ENDIF
+!         
+!     ENDIF 
+     
+
+
+
+
+
+
+
+     
+     
+     IF( lwp ) THEN
+        IF  ( ln_dct_ascii ) THEN
+             !WRITE(numdct_NOOS,'(9e12.4E2)')   zsumclasses( 1)+zsumclasses( 2) ,  zsumclasses( 1), zsumclasses( 2),   &
+             WRITE(numdct_NOOS,'(3F18.3,6e16.8E2)')   zsumclasses( 1)+zsumclasses( 2) ,  zsumclasses( 1), zsumclasses( 2),   &
+                                             zsumclasses( 7)+zsumclasses( 8) ,  zsumclasses( 7), zsumclasses( 8),   &
+                                             zsumclasses( 9)+zsumclasses(10) ,  zsumclasses( 9), zsumclasses(10)
+             CALL FLUSH(numdct_NOOS)
+        ELSE
+             WRITE(numdct_NOOS)   zsumclasses( 1)+zsumclasses( 2) ,  zsumclasses( 1), zsumclasses( 2),   &
+                                  zsumclasses( 7)+zsumclasses( 8) ,  zsumclasses( 7), zsumclasses( 8),   &
+                                  zsumclasses( 9)+zsumclasses(10) ,  zsumclasses( 9), zsumclasses(10)
+             CALL FLUSH(numdct_NOOS) 
+         ENDIF
+     ENDIF
+         
+    IF ( ln_dct_iom_cont ) THEN
+        
+         noos_iom_dummy(:,:,:) = 0.
+         tmp_iom_output(:) = 0.
+         
+         tmp_iom_output(1) = (zsumclasses( 1)+zsumclasses( 2))
+         tmp_iom_output(2) = zsumclasses( 1)
+         tmp_iom_output(3) = zsumclasses( 2)
+         
+         ! Convert to Sv
+         tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
+         tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
+         tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
+         
+         ! limit maximum and minimum values in iom_put
+         if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
+         if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
+         if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
+         if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
+         if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
+         if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
+         
+         ! Set NaN's to Zero         
+         if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
+         if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
+         if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
+         
+         noos_iom_dummy(:,:,1) = tmp_iom_output(1)
+         noos_iom_dummy(:,:,2) = tmp_iom_output(2)
+         noos_iom_dummy(:,:,3) = tmp_iom_output(3)
+         
+         !noos_iom_dummy(:,:,1) = (zsumclasses( 1)+zsumclasses( 2))
+         !noos_iom_dummy(:,:,2) = zsumclasses( 1)
+         !noos_iom_dummy(:,:,3) = zsumclasses( 2)
+         
+         
+         
+         noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_trans'
+         if ( lwp .AND. verbose ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
+         CALL iom_put( noos_var_sect_name,  noos_iom_dummy )
+         noos_iom_dummy(:,:,:) = 0.
+         tmp_iom_output(:) = 0.
+         
+         tmp_iom_output(1) = (zsumclasses( 7)+zsumclasses( 8))
+         tmp_iom_output(2) = zsumclasses( 7)
+         tmp_iom_output(3) = zsumclasses( 8)
+         
+         ! Convert to TJ/s
+         tmp_iom_output(1) = tmp_iom_output(1)*1.E-12
+         tmp_iom_output(2) = tmp_iom_output(2)*1.E-12
+         tmp_iom_output(3) = tmp_iom_output(3)*1.E-12
+         
+         ! limit maximum and minimum values in iom_put
+         if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
+         if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
+         if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
+         if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
+         if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
+         if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
+         
+         ! Set NaN's to Zero         
+         if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
+         if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
+         if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
+         
+         noos_iom_dummy(:,:,1) = tmp_iom_output(1)
+         noos_iom_dummy(:,:,2) = tmp_iom_output(2)
+         noos_iom_dummy(:,:,3) = tmp_iom_output(3)
+         
+         !noos_iom_dummy(:,:,1) = (zsumclasses( 7)+zsumclasses( 8))
+         !noos_iom_dummy(:,:,2) = zsumclasses( 7)
+         !noos_iom_dummy(:,:,3) = zsumclasses( 8)
+         
+         noos_var_sect_name =  "noos_" // trim(noos_sect_name) // '_heat'
+         if ( lwp .AND. verbose ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
+         CALL iom_put(noos_var_sect_name,  noos_iom_dummy )  
+         noos_iom_dummy(:,:,:) = 0.
+         tmp_iom_output(:) = 0.
+         
+         tmp_iom_output(1) = (zsumclasses( 9)+zsumclasses( 10))
+         tmp_iom_output(2) = zsumclasses( 9)
+         tmp_iom_output(3) = zsumclasses( 10)
+         
+         ! Convert to  MT/s
+         tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
+         tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
+         tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
+         
+         
+         ! limit maximum and minimum values in iom_put
+         if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
+         if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
+         if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
+         if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
+         if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
+         if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
+         
+         ! Set NaN's to Zero         
+         if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
+         if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
+         if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
+         
+         noos_iom_dummy(:,:,1) = tmp_iom_output(1)
+         noos_iom_dummy(:,:,2) = tmp_iom_output(2)
+         noos_iom_dummy(:,:,3) = tmp_iom_output(3)
+         
+         !noos_iom_dummy(:,:,1) = (zsumclasses( 9)+zsumclasses( 10))
+         !noos_iom_dummy(:,:,2) = zsumclasses( 9)
+         !noos_iom_dummy(:,:,3) = zsumclasses( 10)
+         
+         noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_salt'
+         if ( lwp .AND. verbose ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
+         CALL iom_put(noos_var_sect_name,  noos_iom_dummy )
+         noos_iom_dummy(:,:,:) = 0.         
+         tmp_iom_output(:) = 0.
+
+
+        DEALLOCATE(noos_iom_dummy)
+     ENDIF
+     
 
      DO jclass=1,MAX(1,sec%nb_class-1)
@@ -1641 +2313 @@
                   
         !write volume transport per class
-        IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN
-           WRITE(numdct_NOOS,'(9e12.4E2)') -(sec%transport( 1,jclass)+sec%transport( 2,jclass)),-sec%transport( 2,jclass),-sec%transport( 1,jclass), &
-                                           -(sec%transport( 7,jclass)+sec%transport( 8,jclass)),-sec%transport( 8,jclass),-sec%transport( 7,jclass), &
-                                           -(sec%transport( 9,jclass)+sec%transport(10,jclass)),-sec%transport(10,jclass),-sec%transport( 9,jclass)
-        ELSE
-           WRITE(numdct_NOOS,'(9e12.4E2)')   sec%transport( 1,jclass)+sec%transport( 2,jclass) , sec%transport( 1,jclass), sec%transport( 2,jclass), &
-                                             sec%transport( 7,jclass)+sec%transport( 8,jclass) , sec%transport( 7,jclass), sec%transport( 8,jclass), &
-                                             sec%transport( 9,jclass)+sec%transport(10,jclass) , sec%transport( 9,jclass), sec%transport(10,jclass)
+        IF( lwp ) THEN
+            
+            IF  ( ln_dct_ascii ) THEN
+                CALL FLUSH(numdct_NOOS)
+
+                !WRITE(numdct_NOOS,'(9e12.4E2)')   sec%transport( 1,jclass)+sec%transport( 2,jclass) , sec%transport( 1,jclass), sec%transport( 2,jclass), &
+                !                                  sec%transport( 7,jclass)+sec%transport( 8,jclass) , sec%transport( 7,jclass), sec%transport( 8,jclass), &
+                !                                  sec%transport( 9,jclass)+sec%transport(10,jclass) , sec%transport( 9,jclass), sec%transport(10,jclass)
+                WRITE(numdct_NOOS,'(3F18.3,6e16.8E2)')   sec%transport( 1,jclass)+sec%transport( 2,jclass) , sec%transport( 1,jclass), sec%transport( 2,jclass), &
+                                                         sec%transport( 7,jclass)+sec%transport( 8,jclass) , sec%transport( 7,jclass), sec%transport( 8,jclass), &
+                                                         sec%transport( 9,jclass)+sec%transport(10,jclass) , sec%transport( 9,jclass), sec%transport(10,jclass)
+            ELSE
+
+                CALL FLUSH(numdct_NOOS)
+                WRITE(numdct_NOOS)   sec%transport( 1,jclass)+sec%transport( 2,jclass) , sec%transport( 1,jclass), sec%transport( 2,jclass), &
+                                     sec%transport( 7,jclass)+sec%transport( 8,jclass) , sec%transport( 7,jclass), sec%transport( 8,jclass), &
+                                     sec%transport( 9,jclass)+sec%transport(10,jclass) , sec%transport( 9,jclass), sec%transport(10,jclass)
+            ENDIF
         ENDIF
 
      ENDDO
+     
+     !IF  ( ln_dct_ascii ) THEN
+        if ( lwp ) CALL FLUSH(numdct_NOOS)
+     !ENDIF
 
      CALL wrk_dealloc(nb_type_class , zsumclasses )  
@@ -1671 +2357 @@
      !!------------------------------------------------------------- 
      !!arguments
-     INTEGER, INTENT(IN)          :: hr          ! hour
+     INTEGER, INTENT(IN)          :: hr          ! hour => effectively kt/12
      TYPE(SECTION), INTENT(INOUT) :: sec         ! section to write   
      INTEGER ,INTENT(IN)          :: ksec        ! section number
@@ -1682 +2368 @@
      !
      REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace 
+     CHARACTER(len=3)      :: noos_sect_name             ! Classname 
+     CHARACTER(len=25)      :: noos_var_sect_name
+     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   noos_iom_dummy
+     INTEGER               :: IERR
+     LOGICAL       ::   verbose     
+     verbose = .false.
+     
      !!------------------------------------------------------------- 
 
+     IF( lwp .AND. verbose ) THEN
+        WRITE(numout,*) " "
+        WRITE(numout,*) "dia_dct_wri_NOOS_h: write transports through section Transect:",ksec-1," at timestep: ", hr
+        WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
+     ENDIF   
+     
      CALL wrk_alloc(nb_type_class , zsumclasses ) 
+     
+     
+     write (noos_sect_name, "(I03)")  ksec
+     
+     ALLOCATE( noos_iom_dummy(jpi,jpj,3),  STAT= ierr )
+        IF( ierr /= 0 )   CALL ctl_stop( 'dia_dct_wri_NOOS_h: failed to allocate noos_iom_dummy array' )
+
+
+
+     
 
      zsumclasses(:)=0._wp
      zslope = sec%slopeSection       
 
+     ! Sum up all classes, to give the total per type (pos vol trans, neg vol trans etc...)
      DO jclass=1,MAX(1,sec%nb_class-1)
         zsumclasses(1:nb_type_class)=zsumclasses(1:nb_type_class)+sec%transport_h(1:nb_type_class,jclass)
      ENDDO
- 
+        
+     
+     ! JT I think changing the sign of output according to the zslope is redundant
+     
      !write volume transport per class
-     IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN
-        z_hr_output(ksec,hr,1)=-(zsumclasses(1)+zsumclasses(2))
-     ELSE
-        z_hr_output(ksec,hr,1)= (zsumclasses(1)+zsumclasses(2))
-     ENDIF
-
+     ! Sum positive and vol trans for all classes in first cell of array
+
+     z_hr_output(ksec,1,1)= (zsumclasses(1)+zsumclasses(2))
+     z_hr_output(ksec,2,1)= zsumclasses(1)
+     z_hr_output(ksec,3,1)= zsumclasses(2)
+
+     ! Sum positive and vol trans for each classes in following cell of array
      DO jclass=1,MAX(1,sec%nb_class-1)
-        IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN
-           z_hr_output(ksec,hr,jclass+1)=-(sec%transport_h(1,jclass)+sec%transport_h(2,jclass))
+        z_hr_output(ksec,1,jclass+1)= (sec%transport_h(1,jclass)+sec%transport_h(2,jclass))
+        z_hr_output(ksec,2,jclass+1)= sec%transport_h(1,jclass)
+        z_hr_output(ksec,3,jclass+1)= sec%transport_h(2,jclass)
+     ENDDO
+
+    
+    IF( lwp )  THEN
+        ! JT IF ( hr .eq. 48._wp ) THEN
+        ! JT    WRITE(numdct_NOOS_h,'(I4,a1,I2,a1,I2,a12,i3,a17,i3)') nyear,'.',nmonth,'.',nday,'   Transect:',ksec-1,'   No. of layers:',sec%nb_class-1
+        ! JT    DO jhr=25,48
+        ! JT       WRITE(numdct_NOOS_h,'(11F12.1)')  z_hr_output(ksec,jhr,1), (z_hr_output(ksec,jhr,jclass+1), jclass=1,MAX(1,10) )
+        ! JT    ENDDO
+        ! JT ENDIF
+
+
+
+        IF ( ln_dct_ascii ) THEN
+            WRITE(numdct_NOOS_h,'(I4,a1,I2,a1,I2,a1,I2,a1,I2,a12,i3,a17,i3)') nyear,'.',nmonth,'.',nday,'.',MOD(hr,24),'.',0,'   Transect:',ksec-1,'   No. of layers:',sec%nb_class-1
+            WRITE(numdct_NOOS_h,'(11F18.3)')  z_hr_output(ksec,1,1), (z_hr_output(ksec,1,jclass+1), jclass=1,MAX(1,10) )
+            WRITE(numdct_NOOS_h,'(11F18.3)')  z_hr_output(ksec,2,1), (z_hr_output(ksec,2,jclass+1), jclass=1,MAX(1,10) )
+            WRITE(numdct_NOOS_h,'(11F18.3)')  z_hr_output(ksec,3,1), (z_hr_output(ksec,3,jclass+1), jclass=1,MAX(1,10) )
+            CALL FLUSH(numdct_NOOS_h)
         ELSE
-           z_hr_output(ksec,hr,jclass+1)= (sec%transport_h(1,jclass)+sec%transport_h(2,jclass))
+            WRITE(numdct_NOOS_h) nyear,nmonth,nday,MOD(hr,24),ksec-1,sec%nb_class-1
+            WRITE(numdct_NOOS_h)  z_hr_output(ksec,1,1), (z_hr_output(ksec,1,jclass+1), jclass=1,MAX(1,10) )
+            WRITE(numdct_NOOS_h)  z_hr_output(ksec,2,1), (z_hr_output(ksec,2,jclass+1), jclass=1,MAX(1,10) )
+            WRITE(numdct_NOOS_h)  z_hr_output(ksec,3,1), (z_hr_output(ksec,3,jclass+1), jclass=1,MAX(1,10) )
+            CALL FLUSH(numdct_NOOS_h)
         ENDIF
-     ENDDO
-
-     IF ( hr .eq. 48._wp ) THEN
-        WRITE(numdct_NOOS_h,'(I4,a1,I2,a1,I2,a12,i3,a17,i3)') nyear,'.',nmonth,'.',nday,'   Transect:',ksec-1,'   No. of layers:',sec%nb_class-1
-        DO jhr=25,48
-           WRITE(numdct_NOOS_h,'(11F12.1)')  z_hr_output(ksec,jhr,1), (z_hr_output(ksec,jhr,jclass+1), jclass=1,MAX(1,10) )
-        ENDDO
-     ENDIF
+
+
+     ENDIF 
+
 
      CALL wrk_dealloc(nb_type_class , zsumclasses )
+     
+     DEALLOCATE(noos_iom_dummy)
+
+
 
   END SUBROUTINE dia_dct_wri_NOOS_h
@@ -1730 +2468 @@
      !! 
      !!        2. Write heat transports in "heat_transport"
-     !!           Unit: Peta W : area * Velocity * T * rhop * Cp * 1.e-15
+     !!           Unit: Peta W : area * Velocity * T * rhau * Cp / 1.e15
      !! 
      !!        3. Write salt transports in "salt_transport"
-     !!           Unit: 10^9 Kg/m^2/s : area * Velocity * S * rhop * 1.e-9 
+     !!           Unit: 10^9 g m^3 / s : area * Velocity * S / 1.e6
      !!
      !!------------------------------------------------------------- 
@@ -1810 +2548 @@
            WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope,  &
                               jclass,classe,zbnd1,zbnd2,&
-                              sec%transport(7,jclass)*1.e-15,sec%transport(8,jclass)*1.e-15, &
-                              ( sec%transport(7,jclass)+sec%transport(8,jclass) )*1.e-15
+                              sec%transport(7,jclass)*1000._wp*rcp/1.e15,sec%transport(8,jclass)*1000._wp*rcp/1.e15, &
+                              ( sec%transport(7,jclass)+sec%transport(8,jclass) )*1000._wp*rcp/1.e15
            !write salt transport per class
            WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope,  &
                               jclass,classe,zbnd1,zbnd2,&
-                              sec%transport(9,jclass)*1.e-9,sec%transport(10,jclass)*1.e-9,&
-                              (sec%transport(9,jclass)+sec%transport(10,jclass))*1.e-9
+                              sec%transport(9,jclass)*1000._wp/1.e9,sec%transport(10,jclass)*1000._wp/1.e9,&
+                              (sec%transport(9,jclass)+sec%transport(10,jclass))*1000._wp/1.e9
         ENDIF
 
@@ -1835 +2573 @@
         WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope, &
                            jclass,"total",zbnd1,zbnd2,&
-                           zsumclasses(7)*1.e-15,zsumclasses(8)*1.e-15,&
-                           (zsumclasses(7)+zsumclasses(8) )*1.e-15
+                           zsumclasses(7)* 1000._wp*rcp/1.e15,zsumclasses(8)* 1000._wp*rcp/1.e15,&
+                           (zsumclasses(7)+zsumclasses(8) )* 1000._wp*rcp/1.e15
         !write total salt transport
         WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope, &
                            jclass,"total",zbnd1,zbnd2,&
-                           zsumclasses(9)*1.e-9,zsumclasses(10)*1.e-9,&
-                           (zsumclasses(9)+zsumclasses(10))*1.e-9
+                           zsumclasses(9)*1000._wp/1.e9,zsumclasses(10)*1000._wp/1.e9,&
+                           (zsumclasses(9)+zsumclasses(10))*1000._wp/1.e9
      ENDIF
 
@@ -1878 +2616 @@
   !!    |    I          |    I+1           |    Z=temperature/salinity/density at U-poinT
   !!    |               |                  |
-  !!  ----------------------------------------  1. Veritcal interpolation: compute zbis
+  !!  ----------------------------------------  1. Veritcale interpolation: compute zbis
   !!    |               |                  |       interpolation between ptab(I,J,K) and ptab(I,J,K+1)
   !!    |               |                  |       zbis = 

branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90

@@ -11 +11 @@
    USE iom             ! I/0 library
    USE wrk_nemo        ! working arrays
+   USE diaregmean
 #if defined key_fabm
    USE trc, ONLY: trn

branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90

@@ -48 +48 @@
    USE dia25h          ! 25h Mean output
    USE diaopfoam       ! Diaopfoam output
+   USE diaregmean      ! regionalmean
+   USE diapea          ! pea
    USE iom
    USE ioipsl
@@ -406 +408 @@
       IF (ln_diaopfoam) THEN
          CALL dia_diaopfoam
+      ENDIF
+      if ( ln_pea ) THEN
+         CALL dia_pea( kt )
+      ENDIF
+      IF (ln_diaregmean) THEN
+         CALL dia_regmean( kt )
       ENDIF
       !

branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90

@@ -56 +56 @@
    PUBLIC iom_getatt, iom_use, iom_context_finalize
 
+   INTEGER , PUBLIC ::   n_regions_output
+
    PRIVATE iom_rp0d, iom_rp1d, iom_rp2d, iom_rp3d
    PRIVATE iom_g0d, iom_g1d, iom_g2d, iom_g3d, iom_get_123d
@@ -106 +108 @@
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z_bnds
       !!----------------------------------------------------------------------
+      
+      
+      
+      
+      
+      REAL(wp),  ALLOCATABLE,   DIMENSION(:,:) ::   tmpregion !: temporary region_mask
+      INTEGER, DIMENSION(3) ::   zdimsz   ! number of elements in each of the 3 dimensions (i.e., lon, lat, no of masks, 297,  375,  4) for an array
+      INTEGER               ::   zndims   ! number of dimensions in an array (i.e. 3, )
+      INTEGER :: inum, nmasks,ierr,maskno,idmaskvar,tmpint
+      REAL(wp), ALLOCATABLE,   DIMENSION(:,:,:)  ::   tmp_region_mask_real   ! tempory region_mask of reals
+      
+      LOGICAL ::   ln_diaregmean  ! region mean calculation
+   
+    
+      INTEGER ::   ios                  ! Local integer output status for namelist read
+      LOGICAL :: ln_diaregmean_ascii  ! region mean calculation ascii output
+      LOGICAL :: ln_diaregmean_bin  ! region mean calculation binary output
+      LOGICAL :: ln_diaregmean_nc  ! region mean calculation netcdf output
+      LOGICAL :: ln_diaregmean_karamld  ! region mean calculation including kara mld terms
+      LOGICAL :: ln_diaregmean_pea  ! region mean calculation including pea terms
+      LOGICAL :: ln_diaregmean_diaar5  ! region mean calculation including AR5 SLR terms
+      LOGICAL :: ln_diaregmean_diasbc  ! region mean calculation including Surface BC
+    
+#if defined key_fabm
+      LOGICAL :: ln_diaregmean_bgc  ! region mean calculation including BGC
+#endif
+      ! Read the number region mask to work out how many regions are needed.
+      
+#if defined key_fabm
+      NAMELIST/nam_diaregmean/ ln_diaregmean,ln_diaregmean_ascii,ln_diaregmean_bin,ln_diaregmean_nc,&
+        & ln_diaregmean_karamld, ln_diaregmean_pea,ln_diaregmean_diaar5,ln_diaregmean_diasbc,ln_diaregmean_bgc
+#else
+      NAMELIST/nam_diaregmean/ ln_diaregmean,ln_diaregmean_ascii,ln_diaregmean_bin,ln_diaregmean_nc,&
+        & ln_diaregmean_karamld, ln_diaregmean_pea,ln_diaregmean_diaar5,ln_diaregmean_diasbc
+#endif
+      
+      ! read in Namelist. 
+      !!----------------------------------------------------------------------
+      !
+      REWIND ( numnam_ref )              ! Read Namelist nam_diatmb in referdiatmbence namelist : TMB diagnostics
+      READ   ( numnam_ref, nam_diaregmean, IOSTAT=ios, ERR= 901 )
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaregmean in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist nam_diatmb in configuration namelist  TMB diagnostics
+      READ  ( numnam_cfg, nam_diaregmean, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaregmean in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, nam_diaregmean )
+
+      IF (ln_diaregmean) THEN
+      
+        ! Open region mask for region means, and retrieve the size of the mask (number of levels)          
+          CALL iom_open ( 'region_mask.nc', inum )
+          idmaskvar = iom_varid( inum, 'mask', kdimsz=zdimsz, kndims=zndims, ldstop = .FALSE.)          
+          nmasks = zdimsz(3)
+          
+          ! read in the region mask (which contains floating point numbers) into a temporary array of reals.
+          ALLOCATE( tmp_region_mask_real(jpi,jpj,nmasks),  STAT= ierr )
+          IF( ierr /= 0 )   CALL ctl_stop( 'dia_regmean_init: failed to allocate tmp_region_mask_real array' )
+          
+          ! Use jpdom_unknown to read in a n layer mask.
+          tmp_region_mask_real(:,:,:) = 0
+          CALL iom_get( inum, jpdom_unknown, 'mask', tmp_region_mask_real(1:nlci,1:nlcj,1:nmasks),   &
+              &          kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,nmasks /) )
+          
+          CALL iom_close( inum )
+          !Convert the region mask of reals into one of integers. 
+          
+          
+          n_regions_output = 0
+          DO maskno = 1,nmasks
+              tmpint = maxval(int(tmp_region_mask_real(:,:,maskno)))
+              CALL mpp_max( tmpint )
+              n_regions_output = n_regions_output + (tmpint + 1)
+          END DO
+      
+          
+        
+      ELSE
+        n_regions_output = 1
+      ENDIF
+      
+      
+      
 #if ! defined key_xios2
       ALLOCATE( z_bnds(jpk,2) )
@@ -227 +312 @@
       CALL iom_set_axis_attr( "iax_20C", (/ REAL(20,wp) /) )
       CALL iom_set_axis_attr( "iax_28C", (/ REAL(28,wp) /) )
+      
+      
+      
+      CALL iom_set_axis_attr( "region", (/ (REAL(ji,wp), ji=1,n_regions_output) /) )
+
+      CALL iom_set_axis_attr( "noos", (/ (REAL(ji,wp), ji=1,3) /) )
+
       
       ! automatic definitions of some of the xml attributs

branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/OBS/obs_level_search.h90

@@ -21 +21 @@
       !!        !  2006-10  (A. Weaver) Cleanup
       !!        !  2008-10  (K. Mogensen) Remove assumptions on grid.
+      !!        !  2019-07  (R. Renshaw) Don't exceed model bottom.
       !!----------------------------------------------------------------------
 
@@ -46 +47 @@
             IF ( pgrddep(jk) >= pobsdep(ji) ) EXIT depk
          END DO depk
-         kobsk(ji) = jk
+         kobsk(ji) = MIN( jk, kgrd )
       END DO
 

branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_z1d.h90

@@ -26 +26 @@
       !!      ! 06-10 (A. Weaver) Cleanup
       !!      ! 07-01 (K. Mogensen) Use profile rather than single level
+      !!      ! 19-07 (R. Renshaw) Avoid division by zero
       !!---------------------------------------------------------------------
 
@@ -62 +63 @@
          z1dm = ( pdep(kkco(jdep)) - pobsdep(jdep)      )
          z1dp = ( pobsdep(jdep)    - pdep(kkco(jdep)-1) )
+!        Where ob is below model bottom, use model bottom rather than extrapolate
+         IF ( pdep(kkco(jdep)) <= pobsdep(jdep) ) z1dm = 0.0_wp
+!        Where lower level is missing, use higher level
          IF ( pobsmask(kkco(jdep)) == 0.0_wp ) z1dp = 0.0_wp
 
          zsum = z1dm + z1dp
+         ! if pobsmask==0 and model level depth==observed depth, we get zsum=0
+         IF ( zsum > 0.0_wp ) THEN
          
          IF ( k1dint == 0 ) THEN
@@ -88 +94 @@
 
          ENDIF
+
+         ELSE  ! take value directly from the higher model level
+            pobs(jdep)  = pobsk(kkco(jdep)-1)
+         ENDIF
+
       END DO
 

branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90

@@ -88 +88 @@
    USE stopts
    USE diatmb          ! Top,middle,bottom output
+   USE diaregmean      ! regional means output
+   USE diapea          ! potential energy anomaly output
    USE dia25h          ! 25h mean output
    USE diaopfoam       ! FOAM operational output
@@ -506 +508 @@
       IF(lwp) WRITE(numout,*) 'Euler time step switch is ', neuler
                             CALL dia_tmb_init  ! TMB outputs
+                            CALL dia_regmean_init  ! TMB outputs
+                            CALL dia_pea_init  ! TMB outputs
                             CALL dia_25h_init  ! 25h mean  outputs
                             CALL dia_diaopfoam_init  ! FOAM operational output
@@ -643 +647 @@
       IF( numdct_heat     /= -1 )   CLOSE( numdct_heat     )   ! heat transports
       IF( numdct_salt     /= -1 )   CLOSE( numdct_salt     )   ! salt transports
+      IF( numdct_NOOS     /= -1 )   CLOSE( numdct_NOOS     )   ! NOOS transports
 
       !

branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/TOP_SRC/trcbdy.F90

@@ -109 +109 @@
             ii = idx%nbi(ib,igrd)
             ij = idx%nbj(ib,igrd)
-            zwgt = idx%nbw(ib,igrd)
+            zwgt = idx%nbw(ib,igrd) * rdttrc(ik) / 86400.d0  ! damping with a timescale of day
             tra(ii,ij,ik,jn) = ( tra(ii,ij,ik,jn) + zwgt * ( ( dta%trc(ib,ik) * dta%rn_fac)  & 
                         &  - tra(ii,ij,ik,jn) ) ) * tmask(ii,ij,ik)
          END DO
       END DO 
+! tra is overwritten at the boundary so damping doesn't work here - need neumann in addition
+! to duplicate the internal value at the boundary
+      CALL bdy_trc_nmn( jn,  idx, dta, kt )
+
+
       !
       IF( kt .eq. nit000 ) CLOSE( unit = 102 )