Index: branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90	(revision 11635)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/DIA/diaar5.F90	(revision 11639)
@@ -39,4 +39,9 @@
    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,5 +61,7 @@
       !!----------------------------------------------------------------------
       !
-      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,4 +92,9 @@
       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,4 +133,36 @@
       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,10 +191,13 @@
       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,4 +258,5 @@
       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,4 +282,5 @@
       !
       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,6 +302,13 @@
       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,4 +318,5 @@
                   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,4 +325,5 @@
       !
       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')
Index: branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90	(revision 11635)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90	(revision 11639)
@@ -35,4 +35,5 @@
   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,15 +66,22 @@
   !! * 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,5 +142,5 @@
      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,6 +157,9 @@
      !!
      !!---------------------------------------------------------------------
-     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,18 +175,50 @@
 
      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,4 +226,5 @@
            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,28 +239,93 @@
 
      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,5 +340,5 @@
      !!  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,4 +361,6 @@
      REAL(wp), POINTER, DIMENSION(:)    :: zwork !   "  
      REAL(wp), POINTER, DIMENSION(:,:,:):: zsum  !   "
+     LOGICAL       ::   verbose     
+     verbose = .false.
 
      !!---------------------------------------------------------------------    
@@ -266,121 +377,153 @@
      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,4 +562,6 @@
      LOGICAL :: llbon                                       ,&!local logical
                 lldebug                                       !debug the section
+     LOGICAL       ::   verbose     
+     verbose = .false.
      !!-------------------------------------------------------------------------------------
      CALL wrk_alloc( nb_point_max, directemp )
@@ -424,5 +569,20 @@
      !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,25 +593,37 @@
 
         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,5 +641,5 @@
         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,5 +655,5 @@
 
             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,5 +672,5 @@
            !read points'coordinates and directions 
            !--------------------------------------
-           IF ( ln_NOOS ) THEN
+           IF ( ln_NOOS .AND. verbose ) THEN
               WRITE(numout,*) 'Coords and direction read'
            ENDIF
@@ -558,10 +730,10 @@
            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,4 +774,6 @@
  
      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,6 +877,6 @@
      !!              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,15 +887,15 @@
     
      !! * 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,4 +904,18 @@
      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,9 +939,4 @@
         !                                                      
         !----------------------------------------------------------------------------------------------------
-        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,31 +946,62 @@
         !--------------------------------------!
         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,57 +1014,57 @@
            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,24 +1094,13 @@
    
               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,4 +1160,18 @@
 
         !----------------------------------------------------------------------------------------------------
+        !----------------------------------------------------------------------------------------------------
+        !----------------------------------------------------------------------------------------------------
+        !
+        !
+        ! ! ! ! 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,8 +1194,4 @@
         !                                                      
         !----------------------------------------------------------------------------------------------------
-        isgnu = 1
-        IF( sec%slopeSection .GT. 0 ) THEN  ; isgnv = -1 
-        ELSE                                ; isgnv =  1
-        ENDIF
 
         IF( ld_debug )write(numout,*)"isgnu isgnv ",isgnu,isgnv
@@ -988,31 +1201,62 @@
         !--------------------------------------!
         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,8 +1275,8 @@
            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,215 +1395,5 @@
      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,14 +1450,229 @@
               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,15 +1721,19 @@
                     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,12 +1825,38 @@
      !
      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,14 +1866,407 @@
      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,15 +2313,29 @@
                   
         !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,5 +2357,5 @@
      !!------------------------------------------------------------- 
      !!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,38 +2368,90 @@
      !
      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,8 +2468,8 @@
      !! 
      !!        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,11 +2548,11 @@
            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,11 +2573,11 @@
         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,5 +2616,5 @@
   !!    |    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 = 
Index: branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/DIA/diapea.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/DIA/diapea.F90	(revision 11639)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/DIA/diapea.F90	(revision 11639)
@@ -0,0 +1,299 @@
+MODULE diapea
+   !!======================================================================
+   !!                       ***  MODULE  diapea  ***
+   !! Potential Energy Anomaly 
+   !!======================================================================
+   !! History :  3.6  !  12/2016  (J Tinker)  Original code
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain
+   USE in_out_manager  ! I/O units
+   USE iom             ! I/0 library
+   USE wrk_nemo        ! working arrays
+   USE eosbn2          ! Equation of state - in situ and potential density
+   USE phycst          ! physical constant
+
+   IMPLICIT NONE
+   PRIVATE 
+   
+   PUBLIC   dia_pea_init            ! routine called by nemogcm.F90
+   PUBLIC   dia_pea                 ! routine called by diawri.F90
+   REAL(wp), PUBLIC, SAVE, ALLOCATABLE,   DIMENSION(:,:)  ::   pea,peat,peas   
+   REAL(wp), SAVE, ALLOCATABLE,   DIMENSION(:,:,:)  ::   wgt_co_mat   ! Weighting array for proportion of grid shallower than cut off depth
+   REAL(wp), SAVE, ALLOCATABLE,   DIMENSION(:,:) ::   t_zmean, s_zmean  !Depth mean temperature and salinity: 2d fields
+   REAL(wp), SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   t_zmean_mat, s_zmean_mat  !Depth mean temperature and salinity: 3d fields
+   REAL(wp) ::   zcutoff
+   LOGICAL , PUBLIC ::   ln_pea  ! region mean calculation
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.6 , NEMO Consortium (2014)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dia_pea_init
+      ! Local variables
+      INTEGER :: ji,jj,jk  ! Dummy loop indices
+      REAL(wp) :: sumz,tmpsumz
+      INTEGER  :: ierr                ! error integer for IOM_get
+      INTEGER ::   ios                  ! Local integer output status for namelist read
+      
+      zcutoff = 200.!200m
+      
+      NAMELIST/nam_pea/ ln_pea
+      
+      
+      !
+      REWIND ( numnam_ref )              ! Read Namelist nam_diatmb in referdiatmbence namelist : TMB diagnostics
+      READ   ( numnam_ref, nam_pea, IOSTAT=ios, ERR= 901 )
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_pea in reference namelist', lwp )
+
+      REWIND( numnam_cfg )              ! Namelist nam_diatmb in configuration namelist  TMB diagnostics
+      READ  ( numnam_cfg, nam_pea, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_pea in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, nam_pea )
+
+      IF(lwp) THEN                   ! Control print
+          WRITE(numout,*)
+          WRITE(numout,*) 'dia_pea_init : Output potential energy anomaly Diagnostics'
+          WRITE(numout,*) '~~~~~~~~~~~~'
+          WRITE(numout,*) 'Namelist nam_pea : set pea output '
+          WRITE(numout,*) 'Switch for pea diagnostics (T) or not (F)  ln_diaregmean  = ', ln_pea
+      ENDIF
+      
+      
+   ALLOCATE( pea(jpi,jpj),  STAT= ierr )
+   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate pea array' )
+   
+   ALLOCATE( peat(jpi,jpj),  STAT= ierr )
+   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate peat array' )
+   
+   ALLOCATE( peas(jpi,jpj),  STAT= ierr )
+   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate peas array' )   
+      
+   ALLOCATE( t_zmean_mat(jpi,jpj,jpk),  STAT= ierr )
+   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate t_zmean_mat array' )
+   
+   ALLOCATE( s_zmean_mat(jpi,jpj,jpk),  STAT= ierr )
+   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate s_zmean_mat array' )
+   
+   ALLOCATE( t_zmean(jpi,jpj),  STAT= ierr )
+   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate t_zmean array' )
+   
+   ALLOCATE( s_zmean(jpi,jpj),  STAT= ierr )
+   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate s_zmean array' )
+   
+   ALLOCATE( wgt_co_mat(jpi,jpj,jpk),  STAT= ierr )
+   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate wgt_co_mat array' )
+   
+   
+   pea(:,:) = 0.
+   peat(:,:) = 0.
+   peas(:,:) = 0.
+    
+   if ( ln_pea ) THEN
+          
+      ! create wgt_co_mat mat, with the proportion of the grid (gdept_0) below cut off (200m)
+      
+      DO jj = 1,jpj
+          DO ji = 1,jpi
+              IF ( tmask(ji,jj,1) == 1.0_wp ) THEN
+                sumz = 0.
+                DO jk = 1,jpk
+                  IF ( tmask(ji,jj,jk) == 1.0_wp ) THEN
+                    tmpsumz = sumz + e3t_n(ji,jj,jk)
+                    IF (sumz .ge. zcutoff) THEN
+                      ! Already too deep
+                      wgt_co_mat(ji,jj,jk) = 0.
+                    ELSE IF (tmpsumz .le. zcutoff) THEN
+                      ! Too shallow
+                      wgt_co_mat(ji,jj,jk) = 1.
+                    ELSE
+                      !proprotion of grid box above cut off depth
+                      wgt_co_mat(ji,jj,jk) = (zcutoff-Sumz)/e3t_n(ji,jj,jk)
+                    END IF
+                    sumz = tmpsumz
+                  endif
+                END DO
+                
+              ELSE
+                !if land, set to 0.
+                DO jk = 1,jpk
+                  wgt_co_mat(ji,jj,jk) = 0.
+                END DO
+                  
+              ENDIF
+          END DO
+       END DO
+   ENDIF
+   
+
+   END SUBROUTINE dia_pea_init
+   
+
+   SUBROUTINE dia_pea(kt)
+   INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+    
+   INTEGER :: ji,jj,jk,ierr  ! Dummy loop indices
+   REAL(wp) :: tmpdenom, tmpnum, maxz
+   !rau0
+   
+   REAL(wp), ALLOCATABLE,   DIMENSION(:,:,:,:) ::   ts_pea_mat,ts_pea_mat_TS_mean,ts_pea_mat_S_mean
+   REAL(wp), ALLOCATABLE,   DIMENSION(:,:,:) ::   tmp_pea_rho,tmp_pea_TS_mean_rho,tmp_pea_S_mean_rho
+   REAL(wp) ::   int_y_pea,int_y_pea_t
+   
+   
+      
+   ALLOCATE( ts_pea_mat(jpi,jpj,jpk,2),  STAT= ierr )
+   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate ts_pea_mat array' )    
+   ALLOCATE( ts_pea_mat_TS_mean(jpi,jpj,jpk,2),  STAT= ierr )
+   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate ts_pea_mat_TS_mean array' )   
+   ALLOCATE( ts_pea_mat_S_mean(jpi,jpj,jpk,2),  STAT= ierr )
+   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate ts_pea_mat_S_mean array' )
+   ALLOCATE( tmp_pea_rho(jpi,jpj,jpk),  STAT= ierr )
+   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate tmp_pea_rho array' )
+   ALLOCATE( tmp_pea_TS_mean_rho(jpi,jpj,jpk),  STAT= ierr )
+   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate tmp_pea_TS_mean_rho array' )
+   ALLOCATE( tmp_pea_S_mean_rho(jpi,jpj,jpk),  STAT= ierr )
+   IF( ierr /= 0 )   CALL ctl_stop( 'dia_pea_init: failed to allocate tmp_pea_S_mean_rho array' )
+
+    pea(:,:)=0.0d0
+    peaT(:,:)=0.0d0
+    peaS(:,:)=0.0d0
+    
+    
+    
+    !calculate the depth mean temperature and salinity of the upper 200m. Save this into a 3d array. Set the value where tmask=0 to be tsn. 
+    
+     DO jj = 1,jpj
+         DO ji = 1,jpi
+            IF ( tmask(ji,jj,1) == 1.0_wp ) THEN  ! if a sea point. 
+            
+             !Depth mean temperature
+             tmpdenom = 0.
+             tmpnum = 0.
+             DO jk = 1,jpk
+                IF ( tmask(ji,jj,jk) == 1.0_wp ) THEN
+                  tmpnum = tmpnum + (wgt_co_mat(ji,jj,jk)*e3t_n(ji,jj,jk)*tsn(ji,jj,jk,jp_tem))
+                  tmpdenom = tmpdenom + (wgt_co_mat(ji,jj,jk)*e3t_n(ji,jj,jk))
+                endif
+             END DO
+             t_zmean(ji,jj) = tmpnum/tmpdenom
+             
+             !Depth mean salinity
+             tmpdenom = 0.
+             tmpnum = 0.
+             DO jk = 1,jpk
+                IF ( tmask(ji,jj,jk) == 1.0_wp ) THEN
+                  tmpnum     = tmpnum + (wgt_co_mat(ji,jj,jk)*e3t_n(ji,jj,jk)*tsn(ji,jj,jk,jp_sal))
+                  tmpdenom = tmpdenom + (wgt_co_mat(ji,jj,jk)*e3t_n(ji,jj,jk))
+                endif
+             END DO
+             s_zmean(ji,jj) = tmpnum/tmpdenom
+             
+             !save into a 3d grid
+             DO jk = 1,jpk
+               IF ( tmask(ji,jj,jk) == 1.0_wp ) THEN
+                 t_zmean_mat(ji,jj,jk) = t_zmean(ji,jj)
+                 s_zmean_mat(ji,jj,jk) = s_zmean(ji,jj)
+               else
+                 t_zmean_mat(ji,jj,jk) = tsn(ji,jj,jk,jp_tem)
+                 s_zmean_mat(ji,jj,jk) = tsn(ji,jj,jk,jp_sal)
+               endif
+             END DO
+            else
+              t_zmean(ji,jj) = tsn(ji,jj,1,jp_tem)
+              s_zmean(ji,jj) = tsn(ji,jj,1,jp_sal)
+              DO jk = 1,jpk
+                t_zmean_mat(ji,jj,jk) = tsn(ji,jj,jk,jp_tem)
+                s_zmean_mat(ji,jj,jk) = tsn(ji,jj,jk,jp_sal)
+              END DO
+            endif
+            
+         END DO
+     END DO
+      
+   !Calculate the density from the depth varying, and depth average temperature and salinity
+   !-----------------------------
+   !-----------------------------
+   
+   ts_pea_mat(:,:,:,:) = tsn(:,:,:,:)
+   
+   ts_pea_mat_TS_mean(:,:,:,1) = t_zmean_mat(:,:,:)
+   ts_pea_mat_TS_mean(:,:,:,2) = s_zmean_mat(:,:,:)
+   
+   ts_pea_mat_S_mean(:,:,:,1) = t_zmean_mat(:,:,:)
+   ts_pea_mat_S_mean(:,:,:,2) = tsn(:,:,:,jp_sal)
+   
+   CALL eos ( ts_pea_mat,         tmp_pea_rho,         gdept_n(:,:,:) )
+   CALL eos ( ts_pea_mat_TS_mean, tmp_pea_TS_mean_rho, gdept_n(:,:,:) )
+   CALL eos ( ts_pea_mat_S_mean,  tmp_pea_S_mean_rho,  gdept_n(:,:,:) )
+   tmp_pea_rho = (tmp_pea_rho * rau0) + rau0
+   tmp_pea_TS_mean_rho = (tmp_pea_TS_mean_rho * rau0) + rau0
+   tmp_pea_S_mean_rho = (tmp_pea_S_mean_rho * rau0) + rau0
+   
+   
+   ! to test the density calculation
+   !CALL iom_put( "tmp_pea_rho" , tmp_pea_rho )                 ! pea
+   !CALL iom_put( "tmp_pea_TS_mean_rho" , tmp_pea_TS_mean_rho )                 ! pea
+   
+   
+   ! Caluclation of the PEA.
+    DO jj = 1,jpj
+      DO ji = 1,jpi
+        pea(ji,jj) = 0.
+        peat(ji,jj) = 0.
+        peas(ji,jj) = 0.
+        maxz = 0.
+        int_y_pea = 0.
+        int_y_pea_t = 0.
+        IF ( tmask(ji,jj,1) == 1.0_wp ) THEN ! for sea points
+        
+          ! the depth integrated calculation is summed up over the depths, and then divided by the depth
+          DO jk = 1,jpk
+            !for each level...
+            
+            IF ( tmask(ji,jj,jk) == 1.0_wp ) THEN ! if above the sea bed...
+              int_y_pea =   -((tmp_pea_TS_mean_rho(ji,jj,jk)) - (tmp_pea_rho(ji,jj,jk)))*9.81*gdept_n(ji,jj,jk)*wgt_co_mat(ji,jj,jk)*e3t_n(ji,jj,jk)
+              int_y_pea_t = -((tmp_pea_S_mean_rho(ji,jj,jk))  - (tmp_pea_rho(ji,jj,jk)))*9.81*gdept_n(ji,jj,jk)*wgt_co_mat(ji,jj,jk)*e3t_n(ji,jj,jk)
+            else
+              int_y_pea = 0.
+              int_y_pea_t = 0.
+            endif
+            
+            ! check that the sum is not NaN. 
+            if ( int_y_pea .ne.  int_y_pea    ) int_y_pea = 0.
+            if ( int_y_pea_t .ne. int_y_pea_t )  int_y_pea_t = 0.
+            !if ( (int_y_pea*int_y_pea    ) .gt. 1.0e6 ) int_y_pea = 0.
+            !if ( (int_y_pea_t*int_y_pea_t) .gt. 1.0e6 ) int_y_pea_t = 0.
+            
+            pea(ji,jj) =  pea(ji,jj) + int_y_pea
+            peat(ji,jj) =  peat(ji,jj) + int_y_pea_t
+            maxz = maxz + (e3t_n(ji,jj,jk)*wgt_co_mat(ji,jj,jk))
+          enddo
+            
+            
+          !divide by the depth
+          pea(ji,jj) = pea(ji,jj)/maxz
+          peat(ji,jj) = peat(ji,jj)/maxz
+          peas(ji,jj) = pea(ji,jj) - peat(ji,jj)
+            
+            
+          else
+            pea(ji,jj) = 0.
+            peat(ji,jj) = 0.
+            peas(ji,jj) = 0.
+          endif
+       enddo
+    enddo
+!    
+     CALL iom_put( "pea" , pea )                 ! pea
+     CALL iom_put( "peat" , peat )               ! pea
+     CALL iom_put( "peas" , peas )               ! pea
+           
+
+    DEALLOCATE(ts_pea_mat,ts_pea_mat_TS_mean,ts_pea_mat_S_mean,tmp_pea_rho,tmp_pea_TS_mean_rho,tmp_pea_S_mean_rho)
+   
+   END SUBROUTINE dia_pea
+
+END MODULE diapea
Index: branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/DIA/diaregmean.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/DIA/diaregmean.F90	(revision 11639)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/DIA/diaregmean.F90	(revision 11639)
@@ -0,0 +1,1111 @@
+MODULE diaregmean 
+   !!======================================================================
+   !!                       ***  MODULE  diaharm  ***
+   !! Timeseries of Regional Means 
+   !!======================================================================
+   !! History :  3.6  !  11/2016  (J Tinker)  Original code
+   !!----------------------------------------------------------------------
+   USE oce             ! ocean dynamics and tracers variables
+   USE dom_oce         ! ocean space and time domain
+   USE in_out_manager  ! I/O units
+   USE iom             ! I/0 library
+   USE wrk_nemo        ! working arrays
+   USE diapea          ! PEA
+   USE zdfmxl          ! MLD
+   USE sbc_oce
+#if defined key_diaar5 
+   USE diaar5
+#endif
+
+
+#if defined key_fabm
+   USE trc
+   USE par_fabm
+#endif
+
+   IMPLICIT NONE
+   PRIVATE
+
+   LOGICAL , PUBLIC ::   ln_diaregmean  ! region mean calculation
+   PUBLIC   dia_regmean_init            ! routine called by nemogcm.F90
+   PUBLIC   dia_regmean                 ! routine called by diawri.F90   
+   PUBLIC   dia_calctmb_region_mean     ! routine called by diatmb.F90
+
+   
+   
+   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_diaar5      ! region mean calculation including AR5 SLR terms
+   LOGICAL :: ln_diaregmean_diasbc      ! region mean calculation including Surface BC
+   LOGICAL :: ln_diaregmean_karamld     ! region mean calculation including kara mld terms
+   LOGICAL :: ln_diaregmean_pea         ! region mean calculation including pea terms
+
+
+   LOGICAL :: ln_diaregmean_bgc         ! region mean calculation including BGC terms
+
+
+   
+   REAL(wp), SAVE, ALLOCATABLE,   DIMENSION(:,:,:)  ::   tmp_region_mask_real   ! tempory region_mask of reals
+   INTEGER,  SAVE, ALLOCATABLE,   DIMENSION(:,:,:)  ::   region_mask            ! region_mask matrix
+   INTEGER                                          ::   nmasks                 ! Number of mask files in region_mask.nc file - 
+   INTEGER,  SAVE, ALLOCATABLE,   DIMENSION(:)      ::   nreg_mat               ! Number of regions in each mask
+   
+   REAL(wp),  ALLOCATABLE,   DIMENSION(:,:,:) ::   tmp_field_mat !: temporary region_mask
+   REAL(wp),  ALLOCATABLE,   DIMENSION(:,:,:) ::   tmp_field_AR5_mat !: temporary region_mask
+   REAL(wp),  ALLOCATABLE,   DIMENSION(:,:,:) ::   tmp_field_SBC_mat !: temporary region_mask
+   INTEGER  ::   tmp_field_cnt                                   ! tmp_field_cnt integer
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.6 , NEMO Consortium (2014)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE dia_regmean_init 
+      !!---------------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_regmean_init  ***
+      !!     
+      !! ** Purpose: Initialization of region mask namelist 
+      !!        
+      !! ** Method : Read namelist
+      !!   History
+      !!   3.6  !  11-16  (J Tinker) Routine to initialize dia_regmean
+      !!---------------------------------------------------------------------------
+      !!
+      INTEGER ::   ios                  ! Local integer output status for namelist read
+      INTEGER  ::  inum                 ! temporary logical unit ! copied from DOM/domzgr.F90
+      INTEGER  ::   ierr                ! error integer for IOM_get
+      INTEGER  ::   idmaskvar           ! output of iom_varid
+      INTEGER  ::   maskno              ! counter for number of masks
+      INTEGER  ::   jj,ji               ! i and j index
+      INTEGER  ::   tmpint              ! temporary integer
+      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, )
+
+
+#if defined key_fabm
+      INTEGER               ::   js,jl,jn, tmp_dummy
+
+      CHARACTER (len=120) ::    tmp_name,tmp_long_name, tmp_unit
+
+      INTEGER               ::   BGC_nlevs,nBGC_output, bgci
+      CHARACTER(len = 10), ALLOCATABLE, DIMENSION(:)       ::   BGC_stat_name(:),BGC_lev_name(:),BGC_output_var(:)
+#endif
+
+      !
+      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
+      
+      
+      ! read in Namelist. 
+      !!----------------------------------------------------------------------
+      !
+      ! Defaults are generally negative, regional means output switched off
+      ln_diaregmean         = .FALSE.  ! don't write regional means
+      ln_diaregmean_bgc     = .FALSE.  ! don't write biogeochemistry regional means
+      ln_diaregmean_karamld = .FALSE.  ! don't write mixed layer depth
+      ln_diaregmean_pea     = .FALSE.  ! don't write potential energy anomaly
+      ln_diaregmean_diaar5  = .FALSE.  ! don't write AR5 SLR terms
+      ln_diaregmean_diasbc  = .FALSE.  ! don't write surface forcing terms
+
+      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(lwp) THEN                   ! Control print
+          WRITE(numout,*)
+          WRITE(numout,*) 'dia_regmean_init : Output regional mean Diagnostics'
+          WRITE(numout,*) '~~~~~~~~~~~~'
+          WRITE(numout,*) 'Namelist nam_regmean : set regmeanoutputs '
+          WRITE(numout,*) 'Switch for regmean diagnostics (T) or not (F)  ln_diaregmean  = ', ln_diaregmean
+          WRITE(numout,*) 'Switch for regmean ascii output (T) or not (F)  ln_diaregmean_ascii  = ', ln_diaregmean_ascii
+          WRITE(numout,*) 'Switch for regmean binary output (T) or not (F)  ln_diaregmean_bin  = ', ln_diaregmean_bin
+          WRITE(numout,*) 'Switch for regmean netcdf output (T) or not (F)  ln_diaregmean_nc  = ', ln_diaregmean_nc
+          WRITE(numout,*) 'Switch for regmean kara mld terms (T) or not (F)  ln_diaregmean_karamld  = ', ln_diaregmean_karamld
+          WRITE(numout,*) 'Switch for regmean PEA terms (T) or not (F)  ln_diaregmean_pea  = ', ln_diaregmean_pea
+          WRITE(numout,*) 'Switch for regmean AR5 SLR terms (T) or not (F)  ln_diaregmean_diaar5  = ', ln_diaregmean_diaar5
+          WRITE(numout,*) 'Switch for regmean Surface forcing terms (T) or not (F)  ln_diaregmean_diasbc  = ', ln_diaregmean_diasbc
+          WRITE(numout,*) 'Switch for regmean BioGeoChemistry terms (T) or not (F)  ln_diaregmean_bgc  = ', ln_diaregmean_bgc
+      ENDIF
+      
+      
+      ALLOCATE( tmp_field_mat(jpi,jpj,19),  STAT= ierr ) !SS/NB/DT/ZA/VA T/S, SSH, MLD, PEA, PEAT, PEAS
+          IF( ierr /= 0 )   CALL ctl_stop( 'tmp_field_mat: failed to allocate tmp_field_mat array' )
+      tmp_field_mat(:,:,:) = 0.
+      tmp_field_cnt = 0
+      
+      IF(ln_diaregmean_diaar5) THEN   
+        ALLOCATE( tmp_field_AR5_mat(jpi,jpj,4),  STAT= ierr ) !SLR terms
+            IF( ierr /= 0 )   CALL ctl_stop( 'tmp_field_AR5_mat: failed to allocate tmp_field_AR5_mat array' )
+        tmp_field_AR5_mat(:,:,:) = 0.
+      ENDIF
+      
+      IF(ln_diaregmean_diasbc) THEN   
+        ALLOCATE( tmp_field_SBC_mat(jpi,jpj,7),  STAT= ierr ) !SBC terms
+            IF( ierr /= 0 )   CALL ctl_stop( 'tmp_field_SBC_mat: failed to allocate tmp_field_SBC_mat array' )
+        tmp_field_SBC_mat(:,:,:) = 0.
+      ENDIF
+
+
+#if defined key_fabm
+      ! as there are so many BGC variables, write out the necessary iodef.xml and field_def.xml entries into ocean.output
+      
+      IF(ln_diaregmean_bgc) THEN   
+            IF(lwp) THEN                   ! Control print
+
+                BGC_nlevs = 5
+                ALLOCATE( BGC_stat_name(6),BGC_lev_name(BGC_nlevs))
+                nBGC_output = 16
+                ALLOCATE( BGC_output_var(nBGC_output))
+
+                BGC_output_var(1) = 'N1_p'
+                BGC_output_var(2) = 'N3_n'
+                BGC_output_var(3) = 'N4_n'
+                BGC_output_var(4) = 'N5_s'
+                BGC_output_var(5) = 'O2_o'
+                BGC_output_var(6) = 'P1_Chl'
+                BGC_output_var(7) = 'P2_Chl'
+                BGC_output_var(8) = 'P3_Chl'
+                BGC_output_var(9) = 'P4_Chl'
+                BGC_output_var(10) = 'P1_c'
+                BGC_output_var(11) = 'P2_c'
+                BGC_output_var(12) = 'P3_c'
+                BGC_output_var(13) = 'P4_c'
+                BGC_output_var(14) = 'Z4_c'
+                BGC_output_var(15) = 'Z5_c'
+                BGC_output_var(16) = 'Z6_c'
+                
+                BGC_stat_name(1) = '_ave'
+                BGC_stat_name(2) = '_tot'
+                BGC_stat_name(3) = '_var'
+                BGC_stat_name(4) = '_cnt'
+                BGC_stat_name(5) = '_reg_id'
+                BGC_stat_name(6) = '_mask_id'
+                BGC_lev_name(1) = 'top'
+                BGC_lev_name(2) = 'bot'
+                BGC_lev_name(3) = 'dif'
+                BGC_lev_name(4) = 'zav'
+                BGC_lev_name(5) = 'vol'
+
+
+                WRITE(numout,*) ''
+                WRITE(numout,*) 'diaregmean BGC field_def.xml entries'
+                WRITE(numout,*) ''
+                
+                
+                DO jn=1,jp_fabm ! State loop
+                    DO js=1,6
+                        DO jl=1,BGC_nlevs
+               
+                            tmp_name=TRIM( TRIM("reg_")//TRIM(BGC_lev_name(jl))//TRIM("_")//TRIM(ctrcnm(jn))// TRIM(BGC_stat_name(js)) )
+
+                            tmp_long_name  = TRIM(ctrcln(jn))
+                            tmp_unit       = TRIM(ctrcun(jn))
+                            
+                            ! Where using volume integrated values, change units... 
+
+                            IF ((jl .EQ. 5) .AND. (js .EQ. 2)) then
+                                SELECT CASE (trim(tmp_unit))
+                                    CASE ('mg C/m^3') ;     tmp_unit = 'Mg C (T C)' !'mg C/m^3'
+                                    CASE ('mg/m^3') ;       tmp_unit = 'Mg (T)'     !'mg/m^3'
+                                    CASE ('mmol C/m^3') ;   tmp_unit = 'Mmol C'     !'mmol C/m^3'
+                                    CASE ('mmol N/m^3') ;   tmp_unit = 'Mmol N'     !'mmol N/m^3'
+                                    CASE ('mmol O_2/m^3') ; tmp_unit = 'Mmol O'     !'mmol O_2/m^3'
+                                    CASE ('mmol P/m^3') ;   tmp_unit = 'Mmol P'     !'mmol P/m^3'
+                                    CASE ('mmol Si/m^3') ;  tmp_unit = 'Mmol S'     !'mmol Si/m^3'
+                                    CASE ('umol/kg') ;      tmp_unit = 'Mmol'       !'umol/kg'   =  mmol/m^3
+                                   ! CASE ('1/m') ;      cycle
+                                    CASE DEFAULT
+                                           tmp_unit = TRIM(TRIM(tmp_unit)//TRIM('x 1e9 m^3'))
+                                END SELECT
+                            ENDIF
+
+                            WRITE(numout,*) TRIM(TRIM('<field id="')//TRIM(tmp_name)//TRIM('"         long_name="')// &
+         &                        TRIM(BGC_lev_name(jl))//TRIM('_')//TRIM(tmp_long_name)//TRIM(BGC_stat_name(js))// &
+         &                        TRIM('"     unit="'//TRIM(tmp_unit)  //'" />'))
+
+                        END DO
+                    END DO
+                END DO
+        
+                WRITE(numout,*) ''
+                WRITE(numout,*) 'diaregmean BGC iodef.xml entries'
+                WRITE(numout,*) ''
+                DO js=1,6
+
+                    DO jn=1,jp_fabm ! State loop
+                        
+                        DO bgci=1,nBGC_output!
+                            if (trim(ctrcnm(jn)) == TRIM(BGC_output_var(bgci))) CYCLE
+                        ENDDO
+                        DO jl=1,BGC_nlevs
+                            ! only print out area averages for ss, nb, diff, and depth averaged, and total values for volume integrated
+                            IF ((jl .EQ. 5) .AND. (js .NE. 2)) CYCLE ! cycle if vol, and not tot.
+                            IF ((jl .NE. 5) .AND. (js .NE. 1)) CYCLE ! cycle if other levels, and not ave. 
+               
+                            tmp_name=TRIM(TRIM("reg_")//TRIM(BGC_lev_name(jl))//TRIM("_")//TRIM(ctrcnm(jn))// TRIM(BGC_stat_name(js))) 
+                            tmp_long_name  = TRIM(ctrcln(jn))
+
+                            WRITE(numout,*) TRIM(TRIM('<field field_ref="')//TRIM(tmp_name)//TRIM('"/>'))
+
+                        END DO !level
+                    END DO ! State loop
+                END DO !statistic
+                WRITE(numout,*) ''
+                DEALLOCATE( BGC_stat_name,BGC_lev_name)
+
+            ENDIF   ! Control print
+
+      ENDIF !ln_diaregmean_bgc
+      
+#endif
+
+      
+      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. 
+          
+          ALLOCATE( region_mask(jpi,jpj,nmasks),  STAT= ierr )
+          IF( ierr /= 0 )   CALL ctl_stop( 'dia_regmean_init: failed to allocate region_mask array' )
+          region_mask(:,:,:) = 0
+          region_mask = int(tmp_region_mask_real(:,:,:))
+          DEALLOCATE( tmp_region_mask_real)
+          
+          
+          ALLOCATE( nreg_mat(nmasks),  STAT= ierr )
+          IF( ierr /= 0 )   CALL ctl_stop( 'dia_regmean_init: failed to allocate nreg_mat array' )
+
+          ! work out the number of regions in each mask, asssuming land is 0, and the regions are consectively numbered, 
+          ! without missing any number, so the number of regions is the maximum number + 1 (for land). mpp_max across the 
+          ! processors to get the global maxima
+          DO maskno = 1,nmasks
+              tmpint = maxval(region_mask(:,:,maskno))
+              CALL mpp_max( tmpint )
+              nreg_mat(maskno) = tmpint + 1
+          END DO
+        
+          IF(lwp) THEN 
+              ! if writing out as binary and text, open the files. 
+              IF ( ln_diaregmean_bin ) THEN
+                  ! Open binary for region means
+                  CALL ctl_opn( numdct_reg_bin  ,'region_mean_timeseries.dat'  , 'NEW', 'UNFORMATTED', 'SEQUENTIAL', -1, numout,  .TRUE. )
+              ENDIF
+              
+              IF ( ln_diaregmean_ascii ) THEN
+                  ! Open text files for region means
+                  CALL ctl_opn( numdct_reg_txt  ,'region_mean_timeseries.txt'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .TRUE. )
+              ENDIF
+          ENDIF
+     ENDIF
+
+   END SUBROUTINE dia_regmean_init
+
+   SUBROUTINE dia_calctmb_region_mean( pinfield,pouttmb )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_calctmb_region_mean  ***
+      !!                   
+      !! ** Purpose :    Find the Top, Bottom and Top minus Bottom fields of water Column
+      !!            :    and depth average, and volume and mass intergated values. 
+
+      !!
+      !! ** Method  :   
+      !!      use mbathy to find surface, mid and bottom of model levels
+      !!
+      !! History :
+      !!   3.6  !  08-14  (E. O'Dea) Routine based on dia_wri_foam
+      !!----------------------------------------------------------------------
+      !! * Modules used
+
+      ! Routine to map 3d field to top, middle, bottom
+      IMPLICIT NONE
+
+
+      ! Routine arguments
+      REAL(wp), DIMENSION(jpi, jpj, jpk), INTENT(IN   ) :: pinfield    ! Input 3d field and mask
+      REAL(wp), DIMENSION(jpi, jpj, 6  ), INTENT(  OUT) :: pouttmb     ! Output top, bottom and surface minus bed, zav, vol int, mass int
+
+      ! Local variables
+      INTEGER :: ji,jj,jk  ! Dummy loop indices
+
+      ! Local Real
+      REAL(wp)                         ::   zmdi  !  set masked values
+      ! for depth int
+      REAL(wp)                         ::   tmpnumer,tmpnumer_mass,tmpdenom ,z_av_val,vol_int_val
+
+      zmdi=1.e+20 !missing data indicator for masking
+
+      !zmdi=0 !missing data indicator for masking
+
+      ! Calculate top
+      pouttmb(:,:,1) = pinfield(:,:,1)*tmask(:,:,1)  + zmdi*(1.0-tmask(:,:,1))
+
+     ! Calculate middle
+      !DO jj = 1,jpj
+      !    DO ji = 1,jpi
+      !        jk              = max(1,mbathy(ji,jj)/2)
+      !        pouttmb(ji,jj,2) = pinfield(ji,jj,jk)*tmask(ji,jj,jk)  + zmdi*(1.0-tmask(ji,jj,jk))
+      !    END DO
+      !END DO
+
+      ! Calculate bottom, and top minus bottom
+      DO jj = 1,jpj
+          DO ji = 1,jpi
+              IF ( tmask(ji,jj,1) .EQ. 1) THEN ! if land
+               
+                  jk              = max(1,mbathy(ji,jj) - 1)
+                  pouttmb(ji,jj,2) = pinfield(ji,jj,jk)*tmask(ji,jj,jk)  + zmdi*(1.0-tmask(ji,jj,jk))
+
+                  pouttmb(ji,jj,3) = (pouttmb(ji,jj,1) - pouttmb(ji,jj,2))*tmask(ji,jj,1)  + zmdi*(1.0-tmask(ji,jj,1))
+
+                  !Depth and volume integral:
+                  !---------------------------
+                  !Vol int = Concentration * vol of grid box, summed over depth.
+                  !Mass int = Concentration * vol of grid box * density of water, summed over depth.
+                  !Depth Average = Vol int divided by * (vol of grid box summed over depth).
+
+                  tmpnumer = 0.
+                  tmpnumer_mass = 0.
+                  tmpdenom = 0.
+                  DO jk = 1,jpk
+                     tmpnumer = tmpnumer + pinfield(ji,jj,jk)*tmask(ji,jj,jk)*e1t(ji,jj)*e2t(ji,jj)*e3t_n(ji,jj,jk)
+                     tmpnumer_mass = tmpnumer_mass + pinfield(ji,jj,jk)*tmask(ji,jj,jk)*e1t(ji,jj)*e2t(ji,jj)*e3t_n(ji,jj,jk)*rhop(ji,jj,jk)
+                     tmpdenom = tmpdenom +                    tmask(ji,jj,jk)*e1t(ji,jj)*e2t(ji,jj)*e3t_n(ji,jj,jk)
+                  END DO
+                  !z_av_val = tmpnumer/tmpdenom
+                  !vol_int_val = tmpnumer
+                  !mass_int_val = tmpnumer*density
+
+                  pouttmb(ji,jj,4) = tmpnumer/tmpdenom ! depth averaged
+                  pouttmb(ji,jj,5) = tmpnumer          ! Vol integrated 
+                  pouttmb(ji,jj,6) = tmpnumer_mass     ! Mass integrated (for heat and salt calcs)
+              ELSE
+                    pouttmb(ji,jj,1) = zmdi
+                    pouttmb(ji,jj,2) = zmdi
+                    pouttmb(ji,jj,3) = zmdi
+                    pouttmb(ji,jj,4) = zmdi
+                    pouttmb(ji,jj,5) = zmdi
+                    pouttmb(ji,jj,6) = zmdi
+              ENDIF
+          END DO
+      END DO
+
+   END SUBROUTINE dia_calctmb_region_mean
+
+
+   SUBROUTINE dia_regmean( kt ) 
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE dia_regmean  ***
+      !! ** Purpose :   Produce regional mean diagnostics
+      !!
+      !! ** Method  :   calls dia_wri_region_mean to calculate and write the regional means for a number of variables, 
+      !!                (calling dia_calctmb_region_mean where necessary).
+      !!                
+      !!                Closes all text and binary files on last time step
+      !!                
+      !!      
+      !!      
+      !!
+      !! History :
+      !!   3.6  !  11-16  (J. Tinker) 
+      !!         
+      !!--------------------------------------------------------------------
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: tmp1mat   ! temporary array of 1's
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwtmbT    ! temporary T workspace 
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwtmbS    ! temporary S workspace 
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwtmb1    ! temporary density workspace 
+      REAL(wp)                            ::   zmdi       ! set masked values
+      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
+      
+      REAL(wp)                         ::   zdt  ! temporary reals
+      INTEGER                          ::   i_steps, ierr         ! no of timesteps per hour, allocation error index
+      INTEGER                          ::   maskno,jj,ji,jk,jm,nreg ! indices of mask, i and j, and number of regions
+
+#if defined key_fabm
+      INTEGER                          ::  jn ,tmp_dummy     ! set masked values
+      REAL(wp)                         ::  tmp_val           ! tmp value, to allow min and max value clamping (not implemented)
+      INTEGER                          ::  jl 
+      CHARACTER (len=60) ::    tmp_name_bgc_top,tmp_name_bgc_bot,tmp_name_bgc_dif, tmp_name_bgc_zav, tmp_name_bgc_vol
+      CHARACTER (len=60) ::    tmp_output_filename
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: zwtmbBGC    ! temporary BGC workspace 
+
+      LOGICAL       ::   verbose 
+      verbose = .FALSE.
+      tmp_val = 0
+#endif
+      zmdi=1.e+20 !missing data indicator for maskin
+
+      IF (ln_diaregmean) THEN
+        ! If regional mean calculations required by namelist
+        ! -----------------
+        ! identify hourly time steps (not used)
+        zdt = rdt
+        IF( nacc == 1 ) zdt = rdtmin
+
+        IF( MOD( 3600,INT(zdt) ) == 0 ) THEN
+            i_steps = 3600/INT(zdt)
+        ELSE
+            CALL ctl_stop('STOP', 'dia_regmean: timestep must give MOD(3600,rdt) = 0 otherwise no hourly values are possible')
+        ENDIF
+        
+        ! Every time step, add physical, SBC, PEA, MLD terms to create hourly sums.
+        ! Every hour, then hourly sums are divided by the number of timesteps in the hour to make hourly means
+        ! These hourly mean values are then used to caluclate the regional means, and output with IOM. 
+#if defined key_fabm
+        ! BGC values are not averaged up over the hour, but are output as hourly instantaneous values. 
+#endif
+
+        
+        !Extract 2d fields from 3d T and S with dia_calctmb_region_mean
+        CALL wrk_alloc( jpi , jpj, 6 , zwtmbT )
+        CALL wrk_alloc( jpi , jpj, 6 , zwtmbS )
+        CALL wrk_alloc( jpi , jpj, 6 , zwtmb1 )
+            
+        CALL dia_calctmb_region_mean(  tsn(:,:,:,jp_tem),zwtmbT)
+        CALL dia_calctmb_region_mean(  tsn(:,:,:,jp_sal),zwtmbS)
+
+        ! To calc regional mean time series of int vol and mass, run region mean code on array of 1's...
+        !   - then when multplying by volume, gives volume,
+        !   - then when multplying by volume*density, gives mass
+
+        CALL wrk_alloc( jpi , jpj, jpk , tmp1mat )
+        DO jj = 1,jpj
+            DO ji = 1,jpi
+                DO jk = 1,jpk
+                    tmp1mat(ji,jj,jk) = 1
+                END DO
+            END DO
+        END DO
+
+        CALL dia_calctmb_region_mean(  tmp1mat,zwtmb1)
+        CALL wrk_dealloc( jpi , jpj, jpk , tmp1mat )
+ 
+        ! Add 2d fields every time step to the hourly total.
+            
+        tmp_field_mat(:,:,1) = tmp_field_mat(:,:,1) + (zwtmbT(:,:,1)*tmask(:,:,1)) !sst
+        tmp_field_mat(:,:,2) = tmp_field_mat(:,:,2) + (zwtmbT(:,:,2)*tmask(:,:,1)) !nbt
+        tmp_field_mat(:,:,3) = tmp_field_mat(:,:,3) + (zwtmbT(:,:,3)*tmask(:,:,1)) !dft
+
+        tmp_field_mat(:,:,4) = tmp_field_mat(:,:,4) + (zwtmbT(:,:,4)*tmask(:,:,1)) !zat
+        tmp_field_mat(:,:,5) = tmp_field_mat(:,:,5) + (zwtmbT(:,:,5)*tmask(:,:,1)) !vat
+        tmp_field_mat(:,:,6) = tmp_field_mat(:,:,6) + ((zwtmbT(:,:,6)*tmask(:,:,1)*4.2e3))! heat
+
+        tmp_field_mat(:,:,7) = tmp_field_mat(:,:,7) + (zwtmbS(:,:,1)*tmask(:,:,1)) !sss
+        tmp_field_mat(:,:,8) = tmp_field_mat(:,:,8) + (zwtmbS(:,:,2)*tmask(:,:,1)) !nbs
+        tmp_field_mat(:,:,9) = tmp_field_mat(:,:,9) + (zwtmbS(:,:,3)*tmask(:,:,1)) !dfs
+
+        tmp_field_mat(:,:,10) = tmp_field_mat(:,:,10) + (zwtmbS(:,:,4)*tmask(:,:,1)) !zas
+        tmp_field_mat(:,:,11) = tmp_field_mat(:,:,11) + (zwtmbS(:,:,5)*tmask(:,:,1)) !vas
+        tmp_field_mat(:,:,12) = tmp_field_mat(:,:,12) + (zwtmbS(:,:,6)*tmask(:,:,1)) !salt
+
+        tmp_field_mat(:,:,13) = tmp_field_mat(:,:,13) + (zwtmb1(:,:,5)*tmask(:,:,1))!vol
+        tmp_field_mat(:,:,14) = tmp_field_mat(:,:,14) + (zwtmb1(:,:,6)*tmask(:,:,1))!mass
+
+        tmp_field_mat(:,:,15) = tmp_field_mat(:,:,15) + (sshn(:,:)*tmask(:,:,1)) !ssh
+        
+        CALL wrk_dealloc( jpi , jpj, 6 , zwtmbT )
+        CALL wrk_dealloc( jpi , jpj, 6 , zwtmbS )
+        CALL wrk_dealloc( jpi , jpj, 6 , zwtmb1 )
+        
+        IF( ln_diaregmean_karamld  ) THEN
+            tmp_field_mat(:,:,16) = tmp_field_mat(:,:,16) + (hmld_zint(:,:)*tmask(:,:,1)) !mldkara
+        ENDIF
+
+        IF( ln_diaregmean_pea  ) THEN
+            tmp_field_mat(:,:,17) = tmp_field_mat(:,:,17) + (pea(:,:)*tmask(:,:,1))  !pea
+            tmp_field_mat(:,:,18) = tmp_field_mat(:,:,18) + (peat(:,:)*tmask(:,:,1)) !peat
+            tmp_field_mat(:,:,19) = tmp_field_mat(:,:,19) + (peas(:,:)*tmask(:,:,1)) !peas
+        ENDIF
+          
+        IF( ln_diaregmean_diaar5  ) THEN
+            tmp_field_AR5_mat(:,:,1) = tmp_field_AR5_mat(:,:,1) + (sshsteric_mat(:,:)*tmask(:,:,1))
+            tmp_field_AR5_mat(:,:,2) = tmp_field_AR5_mat(:,:,2) + (sshthster_mat(:,:)*tmask(:,:,1))
+            tmp_field_AR5_mat(:,:,3) = tmp_field_AR5_mat(:,:,3) + (sshhlster_mat(:,:)*tmask(:,:,1))
+            tmp_field_AR5_mat(:,:,4) = tmp_field_AR5_mat(:,:,4) + (zbotpres_mat(:,:)*tmask(:,:,1))
+        ENDIF
+        
+        IF( ln_diaregmean_diasbc  ) THEN
+            tmp_field_SBC_mat(:,:,1) = tmp_field_SBC_mat(:,:,1) + ((qsr  + qns)*tmask(:,:,1))
+            tmp_field_SBC_mat(:,:,2) = tmp_field_SBC_mat(:,:,2) + (qsr*tmask(:,:,1))
+            tmp_field_SBC_mat(:,:,3) = tmp_field_SBC_mat(:,:,3) + (qns*tmask(:,:,1))
+            tmp_field_SBC_mat(:,:,4) = tmp_field_SBC_mat(:,:,4) + (emp*tmask(:,:,1))
+            tmp_field_SBC_mat(:,:,5) = tmp_field_SBC_mat(:,:,5) + (wndm*tmask(:,:,1))
+            tmp_field_SBC_mat(:,:,6) = tmp_field_SBC_mat(:,:,6) + (pressnow*tmask(:,:,1))
+            tmp_field_SBC_mat(:,:,7) = tmp_field_SBC_mat(:,:,7) + (rnf*tmask(:,:,1))
+        ENDIF
+
+
+        
+        tmp_field_cnt = tmp_field_cnt + 1
+
+        ! On the hour, calculate hourly means from the hourly total,and process the regional means. 
+        
+        IF ( MOD( kt, i_steps ) == 0 .and. kt .ne. nn_it000 ) THEN
+
+
+            CALL dia_wri_region_mean(kt, "sst" , tmp_field_mat(:,:,1)/real(tmp_field_cnt,wp))
+            CALL dia_wri_region_mean(kt, "nbt" , tmp_field_mat(:,:,2)/real(tmp_field_cnt,wp))
+            CALL dia_wri_region_mean(kt, "dft" , tmp_field_mat(:,:,3)/real(tmp_field_cnt,wp))
+
+            CALL dia_wri_region_mean(kt, "zat" , tmp_field_mat(:,:,4)/real(tmp_field_cnt,wp))
+            CALL dia_wri_region_mean(kt, "vat" , tmp_field_mat(:,:,5)/real(tmp_field_cnt,wp))
+            CALL dia_wri_region_mean(kt, "heat" , tmp_field_mat(:,:,6)/real(tmp_field_cnt,wp)/1e12)
+
+            CALL dia_wri_region_mean(kt, "sss" , tmp_field_mat(:,:,7)/real(tmp_field_cnt,wp))
+            CALL dia_wri_region_mean(kt, "nbs" , tmp_field_mat(:,:,8)/real(tmp_field_cnt,wp))
+            CALL dia_wri_region_mean(kt, "dfs" , tmp_field_mat(:,:,9)/real(tmp_field_cnt,wp))
+
+            CALL dia_wri_region_mean(kt, "zas" , tmp_field_mat(:,:,10)/real(tmp_field_cnt,wp))
+            CALL dia_wri_region_mean(kt, "vas" , tmp_field_mat(:,:,11)/real(tmp_field_cnt,wp))
+            CALL dia_wri_region_mean(kt, "salt" , tmp_field_mat(:,:,12)/real(tmp_field_cnt,wp)/1e12)
+
+            CALL dia_wri_region_mean(kt, "vol" , tmp_field_mat(:,:,13)/real(tmp_field_cnt,wp))
+            CALL dia_wri_region_mean(kt, "mass" , tmp_field_mat(:,:,14)/real(tmp_field_cnt,wp))
+
+            CALL dia_wri_region_mean(kt, "ssh" , tmp_field_mat(:,:,15)/real(tmp_field_cnt,wp))
+
+
+            IF( ln_diaregmean_karamld  ) THEN          
+                CALL dia_wri_region_mean(kt, "mldkara" , tmp_field_mat(:,:,16)/real(tmp_field_cnt,wp)) ! tm
+            ENDIF
+
+            IF( ln_diaregmean_pea  ) THEN          
+                CALL dia_wri_region_mean(kt, "pea"  , tmp_field_mat(:,:,17)/real(tmp_field_cnt,wp))
+                CALL dia_wri_region_mean(kt, "peat" , tmp_field_mat(:,:,18)/real(tmp_field_cnt,wp))
+                CALL dia_wri_region_mean(kt, "peas" , tmp_field_mat(:,:,19)/real(tmp_field_cnt,wp)) ! tmb
+            ENDIF
+
+            tmp_field_mat(:,:,:) = 0.
+
+            IF( ln_diaregmean_diaar5  ) THEN
+
+                CALL dia_wri_region_mean(kt, "ssh_steric" ,      tmp_field_AR5_mat(:,:,1)/real(tmp_field_cnt,wp))
+                CALL dia_wri_region_mean(kt, "ssh_thermosteric", tmp_field_AR5_mat(:,:,2)/real(tmp_field_cnt,wp))
+                CALL dia_wri_region_mean(kt, "ssh_halosteric" ,  tmp_field_AR5_mat(:,:,3)/real(tmp_field_cnt,wp))
+                CALL dia_wri_region_mean(kt, "bot_pres" ,        tmp_field_AR5_mat(:,:,4)/real(tmp_field_cnt,wp))
+                tmp_field_AR5_mat(:,:,:) = 0.
+            ENDIF
+
+            IF( ln_diaregmean_diasbc  ) THEN
+
+                CALL dia_wri_region_mean(kt, "qt"   , tmp_field_SBC_mat(:,:,1)/real(tmp_field_cnt,wp))
+                CALL dia_wri_region_mean(kt, "qsr"  , tmp_field_SBC_mat(:,:,2)/real(tmp_field_cnt,wp))
+                CALL dia_wri_region_mean(kt, "qns"  , tmp_field_SBC_mat(:,:,3)/real(tmp_field_cnt,wp))
+                CALL dia_wri_region_mean(kt, "emp"  , tmp_field_SBC_mat(:,:,4)/real(tmp_field_cnt,wp))
+                CALL dia_wri_region_mean(kt, "wspd" , tmp_field_SBC_mat(:,:,5)/real(tmp_field_cnt,wp))
+                CALL dia_wri_region_mean(kt, "mslp" , tmp_field_SBC_mat(:,:,6)/real(tmp_field_cnt,wp))
+                CALL dia_wri_region_mean(kt, "rnf"  , tmp_field_SBC_mat(:,:,7)/real(tmp_field_cnt,wp))
+                tmp_field_SBC_mat(:,:,:) = 0.
+            ENDIF
+
+#if defined key_fabm
+            !ADD Biogeochemistry
+            
+            IF( ln_diaregmean_bgc  ) THEN !ln_diaregmean_bgc
+                
+                ! Loop through 3d BGC tracers
+                DO jn=1,jp_fabm ! State loop
+
+                    ! get variable name for different levels
+                    tmp_name_bgc_top=TRIM(TRIM("top_")//TRIM(ctrcnm(jn)))
+                    tmp_name_bgc_bot=TRIM(TRIM("bot_")//TRIM(ctrcnm(jn)))
+                    tmp_name_bgc_dif=TRIM(TRIM("dif_")//TRIM(ctrcnm(jn)))
+                    tmp_name_bgc_zav=TRIM(TRIM("zav_")//TRIM(ctrcnm(jn)))
+                    tmp_name_bgc_vol=TRIM(TRIM("vol_")//TRIM(ctrcnm(jn)))
+
+                    ! print out names if verbose
+                    IF(verbose .AND. lwp) THEN 
+                        WRITE(numout,*)
+                        WRITE(numout,*) 'dia_regmean tmp_name_bgc_top : ',TRIM(tmp_name_bgc_top)
+                        WRITE(numout,*) 'dia_regmean tmp_name_bgc_bot : ',TRIM(tmp_name_bgc_bot)
+                        WRITE(numout,*) 'dia_regmean tmp_name_bgc_dif : ',TRIM(tmp_name_bgc_dif)
+                        WRITE(numout,*) 'dia_regmean tmp_name_bgc_zav : ',TRIM(tmp_name_bgc_zav)
+                        WRITE(numout,*) 'dia_regmean tmp_name_bgc_vol : ',TRIM(tmp_name_bgc_vol)
+                        CALL FLUSH(numout) 
+
+                    ENDIF
+
+                    !Allocate working array, and get surface, bed etc fields. 
+                    CALL wrk_alloc( jpi , jpj,  6 , zwtmbBGC )
+                    CALL dia_calctmb_region_mean(  trn(:,:,:,jn),zwtmbBGC )
+
+
+                    !Print out 2d fields to ascii text files to check values if verbose. (24MB per time step, per BGC variable)
+                    IF (verbose) THEN 
+
+                        WRITE (tmp_output_filename, "(A4,I3.3,A1,I6.6,A1,I3.3,A4)") "bgc_",jn,"_",kt,"_",narea,".txt"
+                        WRITE (*,*) tmp_output_filename
+                        OPEN(UNIT=74,FILE=TRIM(tmp_output_filename))
+
+                        DO ji = 1,jpi
+                            DO jj = 1,jpj
+                                WRITE(74,FMT="(I4,I4,F3,F25.5,F25.5,F25.5,F25.5,F25.5)") nimpp+ji, njmpp+jj,tmask(ji,jj,1),&
+                                      & zwtmbBGC(ji,jj,1),zwtmbBGC(ji,jj,2),zwtmbBGC(ji,jj,3),zwtmbBGC(ji,jj,4),zwtmbBGC(ji,jj,5)/1e9
+                            END DO
+                        END DO
+                        CLOSE(74)
+                    ENDIF
+                    
+                    ! Do region means
+                    CALL dia_wri_region_mean(kt, TRIM(tmp_name_bgc_top)   , zwtmbBGC(:,:,1))
+                    CALL dia_wri_region_mean(kt, TRIM(tmp_name_bgc_bot)   , zwtmbBGC(:,:,2))
+                    CALL dia_wri_region_mean(kt, TRIM(tmp_name_bgc_dif)   , zwtmbBGC(:,:,3))
+                    CALL dia_wri_region_mean(kt, TRIM(tmp_name_bgc_zav)   , zwtmbBGC(:,:,4))
+                    CALL dia_wri_region_mean(kt, TRIM(tmp_name_bgc_vol)   , zwtmbBGC(:,:,5)/1e9)
+
+
+                    !Deallocate working array
+                    CALL wrk_dealloc( jpi , jpj,  6 , zwtmbBGC )
+                ENDDO ! State loop
+            ENDIF !ln_diaregmean_bgc
+
+#endif
+              
+            tmp_field_cnt = 0
+  
+        ENDIF ! ( MOD( kt, i_steps ) == 0  .and. kt .ne. nn_it000 )
+        
+        
+        ! If on the last time step, close binary and ascii files. 
+        IF( kt == nitend ) THEN
+            IF(lwp) THEN
+                IF ( ln_diaregmean_bin ) THEN
+                    !Closing binary files for regional mean time series.
+                    CLOSE(numdct_reg_bin)
+                ENDIF
+                IF ( ln_diaregmean_ascii ) THEN
+                    !Closing text files for regional mean time series.
+                    CLOSE(numdct_reg_txt)
+                ENDIF
+
+                DEALLOCATE( region_mask, nreg_mat, tmp_field_mat)
+                IF( ln_diaregmean_diaar5  ) DEALLOCATE( tmp_field_AR5_mat)
+                IF( ln_diaregmean_diasbc  ) DEALLOCATE( tmp_field_SBC_mat)
+            ENDIF
+        ENDIF
+          
+          
+      ELSE
+        CALL ctl_warn('dia_regmean: regmean diagnostic is set to false you should not have seen this')
+      ENDIF
+      
+   END SUBROUTINE dia_regmean
+   
+   
+   SUBROUTINE dia_wri_region_mean(kt, tmp_name,         infield  )
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE dia_tmb  ***
+      !!                   
+      !! ** Purpose :   Calculate and write region mean time series for 2d arrays
+      !!
+      !! ** Method  :   
+      !!      use 
+      !!
+      !! History :
+      !!   ??  !  15/10/2015  (JTinker) Routine taken from old dia_wri_foam
+      !!----------------------------------------------------------------------
+      !! * Modules used
+      !use lib_mpp
+      !use lib_fortr
+      IMPLICIT NONE
+      
+      INTEGER, INTENT(in) ::   kt
+      CHARACTER (len=*) , INTENT(IN   ) ::    tmp_name
+      REAL(wp), DIMENSION(jpi, jpj), INTENT(IN   ) :: infield    ! Input 3d field and mask
+      
+      ! Local variables
+      INTEGER, DIMENSION(jpi, jpj) :: internal_region_mask    ! Input 3d field and mask 
+      REAL(wp), DIMENSION(jpi, jpj) :: internal_infield    ! Internal data field
+      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   zrmet_ave,zrmet_tot,zrmet_var,zrmet_cnt,zrmet_mask_id,zrmet_reg_id  ,zrmet_min,zrmet_max
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zrmet_out
+      REAL(wp), ALLOCATABLE,   DIMENSION(:) ::   ave_mat,tot_mat,num_mat,var_mat,ssq_mat,cnt_mat,reg_id_mat,mask_id_mat ,min_mat,max_mat   !: region_mask
+      
+      REAL(wp)                         ::   zmdi, zrmet_val      ! set masked values
+      INTEGER :: maskno,nreg  ! ocean time-step indexocean time step            
+      INTEGER :: ji,jj,jk,ind,jm ! Dummy loop indices
+      INTEGER :: reg_ind_cnt ! Dummy loop indices
+      
+      INTEGER  ::   ierr      
+      REAL(wp)  :: tmpreal
+      CHARACTER(LEN=180) :: FormatString,nreg_string,tmp_name_iom
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   dummy_zrmet
+      LOGICAL       ::   verbose     
+      verbose = .false.
+
+
+      zmdi=1.e+20 !missing data indicator for maskin
+      
+      !Allocate output arrays for iomput, set to zmdi, and set a region counter = 1
+      ALLOCATE( zrmet_ave(n_regions_output),  STAT= ierr )
+        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_ave array' )
+      ALLOCATE( zrmet_tot(n_regions_output),  STAT= ierr )
+        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_tot array' )
+      ALLOCATE( zrmet_var(n_regions_output),  STAT= ierr )
+        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_var array' )
+      ALLOCATE( zrmet_cnt(n_regions_output),  STAT= ierr )
+        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_cnt array' )
+      ALLOCATE( zrmet_mask_id(n_regions_output),  STAT= ierr )
+        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_mask_id array' )
+      ALLOCATE( zrmet_reg_id(n_regions_output),  STAT= ierr )
+        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_reg_id array' )
+
+
+      ALLOCATE( zrmet_min(n_regions_output),  STAT= ierr )
+        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_min array' )
+      ALLOCATE( zrmet_max(n_regions_output),  STAT= ierr )
+        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_max array' )
+      
+      ALLOCATE( zrmet_out(jpi,jpj,n_regions_output),  STAT= ierr )
+        IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate zrmet_reg_id array' )
+
+  
+      
+        IF(lwp .AND. verbose) THEN 
+              WRITE(numout,*)
+              WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//';'
+              WRITE(numout,*)
+        ENDIF
+        
+      DO ji = 1,jpi
+          DO jj = 1,jpj
+            internal_infield(ji,jj) = infield(ji,jj)
+          END DO
+      END DO
+        
+      ! Check for NANS # JT   03/09/2018
+      DO ji = 1,jpi
+          DO jj = 1,jpj
+              IF ( tmask(ji,jj,1) == 1.0_wp ) THEN
+                  IF ( internal_infield(ji,jj) .ne. internal_infield(ji,jj) ) THEN
+                      WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//' Nan at (kt,i,j): ',kt,ji - (-jpizoom+1-nimpp+1),jj - (-jpjzoom+1-njmpp+1)
+                      internal_infield(ji,jj) = 0.
+                  ENDIF
+              ELSE                
+                  IF ( internal_infield(ji,jj) .ne. internal_infield(ji,jj) ) THEN
+                      WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//' Masked Nan at (kt,i,j): ',kt,ji - (-jpizoom+1-nimpp+1),jj - (-jpjzoom+1-njmpp+1)
+                      internal_infield(ji,jj) = 0.
+                  ENDIF
+              ENDIF
+          END DO
+      END DO
+      
+      
+      zrmet_ave(:) = zmdi
+      zrmet_tot(:) = zmdi
+      zrmet_var(:) = zmdi
+      zrmet_cnt(:) = zmdi
+      zrmet_mask_id(:) = zmdi
+      zrmet_reg_id(:) = zmdi
+      
+      zrmet_min(:) = zmdi
+      zrmet_max(:) = zmdi
+      reg_ind_cnt = 1
+      
+      
+      ! loop though the masks
+      DO maskno = 1,nmasks
+          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; begin mask loops: ',maskno
+          
+          
+          ! For each mask, get the number of regions (nreg), and a local copy of the region. 
+          nreg = nreg_mat(maskno)
+          internal_region_mask = region_mask(:,:,maskno)
+          
+          ! allocate temporary stat arrays, and set to zero
+          ALLOCATE( ave_mat(nreg),  STAT= ierr )
+          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate ave_mat array' )
+          ALLOCATE( tot_mat(nreg),      STAT= ierr )
+          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate tot_mat array' )
+          ALLOCATE( num_mat(nreg),  STAT= ierr )
+          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate num_mat array' )
+          ALLOCATE( var_mat(nreg),  STAT= ierr )
+          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate var_mat array' )
+          ALLOCATE( ssq_mat(nreg),  STAT= ierr )
+          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate ssq_mat array' )
+          ALLOCATE( cnt_mat(nreg),  STAT= ierr )
+          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate cnt_mat array' )
+
+          ALLOCATE( min_mat(nreg),  STAT= ierr )
+          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate min_mat array' )
+          ALLOCATE( max_mat(nreg),  STAT= ierr )
+          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate max_mat array' )
+
+          ALLOCATE( reg_id_mat(nreg),  STAT= ierr )
+          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate reg_id_mat array' )
+          ALLOCATE( mask_id_mat(nreg),  STAT= ierr )
+          IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate mask_id_mat array' )
+
+
+          
+          ave_mat(:) = 0.
+          tot_mat(:) = 0.
+          num_mat(:) = 0.
+          var_mat(:) = 0.
+          cnt_mat(:) = 0.
+          ssq_mat(:) = 0.
+
+          min_mat(:) = zmdi
+          max_mat(:) = -zmdi
+          reg_id_mat(:) = 0.
+          mask_id_mat(:) = 0.
+          
+          ! loop though the array. for each sea grid box where tmask == 1), 
+          ! read which region the grid box is in, add the value of the gridbox (and its square) 
+          ! to the total for that region, and then increment the counter for that region.
+          !CALL cpu_time(start_reg_mean_loop)
+          !WRITE(numout,*) kt,start_reg_mean_loop
+          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; begin spatial loops: '
+          DO ji = nldi,nlei
+              DO jj = nldj,nlej
+                    IF ( tmask(ji,jj,1) == 1.0_wp ) THEN
+                        ind = internal_region_mask(ji,jj)+1
+                        tot_mat(ind) = tot_mat(ind) + (internal_infield(ji,jj))
+                        ssq_mat(ind) = ssq_mat(ind) + ( internal_infield(ji,jj) *  internal_infield(ji,jj))
+                        cnt_mat(ind) = cnt_mat(ind) + 1.
+
+                        min_mat(ind) = min(min_mat(ind),internal_infield(ji,jj))
+                        max_mat(ind) = max(max_mat(ind),internal_infield(ji,jj))
+                    ENDIF
+              END DO
+          END DO
+          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finish spatial loops: '
+          ! sum the totals, the counts, and the squares across the processors          
+          CALL mpp_sum( tot_mat,nreg )
+          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_sum 1'
+          CALL mpp_sum( ssq_mat,nreg )
+          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_sum 2'
+          CALL mpp_sum( cnt_mat,nreg )
+          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_sum 2'
+
+          CALL mpp_min( min_mat,nreg )
+          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_min'
+          CALL mpp_max( max_mat,nreg )
+          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_wri_region_mean : '//tmp_name//'; finished mpp_max'
+          
+          
+          !calculate the mean and variance from the total, sum of squares and the count. 
+          
+          ave_mat = tot_mat(:)/cnt_mat(:)
+          var_mat = ssq_mat(:)/cnt_mat(:) - (ave_mat(:)*ave_mat(:))
+          
+          
+          !mask array of mask and region number. 
+          DO jj = 1,nreg
+              reg_id_mat(jj) = real(jj-1)
+              mask_id_mat(jj) = real(maskno)
+          END DO
+          
+          
+          !write text and binary, and note region statistics for current mask for later iom_put
+          IF( lwp ) THEN 
+          
+              !Write out ascii and binary if requred
+              IF ( ln_diaregmean_bin ) THEN
+                  !Writing out regional mean time series to binary files
+                  WRITE(numdct_reg_bin) tmp_name,kt,maskno,n_regions_output
+                  WRITE(numdct_reg_bin) ave_mat
+                  WRITE(numdct_reg_bin) tot_mat
+                  WRITE(numdct_reg_bin) var_mat
+                  WRITE(numdct_reg_bin) ssq_mat
+                  WRITE(numdct_reg_bin) cnt_mat
+                  WRITE(numdct_reg_bin) min_mat
+                  WRITE(numdct_reg_bin) max_mat
+              ENDIF
+              
+              IF ( ln_diaregmean_ascii  ) THEN
+                  !Writing out regional mean time series to text files
+
+                  WRITE(nreg_string, "(I5)") nreg
+                  FormatString = "(A30,"//trim(nreg_string)//"F25.3)"
+                  WRITE(numdct_reg_txt, FMT="(A30,I6,I6)") tmp_name,kt,maskno            
+                  WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"ave_mat:", ave_mat
+                  WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"tot_mat:", tot_mat
+                  WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"var_mat:", var_mat
+                  WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"ssq_mat:", ssq_mat
+                  WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"cnt_mat:", cnt_mat
+                  WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"min_mat:", min_mat
+                  WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"max_mat:", max_mat
+                  WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"reg_mat:", reg_id_mat
+                  WRITE(numdct_reg_txt, FMT=trim(FormatString)) trim(tmp_name)//" "//"msk_mat:", mask_id_mat
+
+              ENDIF
+              
+              DO jm = 1,nreg
+                  zrmet_ave(    reg_ind_cnt) =     ave_mat(jm)
+                  zrmet_tot(    reg_ind_cnt) =     tot_mat(jm)
+                  zrmet_var(    reg_ind_cnt) =     var_mat(jm)
+                  zrmet_cnt(    reg_ind_cnt) =     cnt_mat(jm)
+                  zrmet_min(    reg_ind_cnt) =     min_mat(jm)
+                  zrmet_max(    reg_ind_cnt) =     max_mat(jm)
+                  zrmet_reg_id( reg_ind_cnt) =  reg_id_mat(jm)
+                  zrmet_mask_id(reg_ind_cnt) = mask_id_mat(jm)
+                
+                  reg_ind_cnt = reg_ind_cnt + 1 
+              END DO
+          
+          ENDIF
+        
+          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean about to deallocated arrays for ',kt,maskno
+          DEALLOCATE(ave_mat,tot_mat,num_mat,var_mat,ssq_mat,cnt_mat,min_mat,max_mat,reg_id_mat,mask_id_mat)
+
+          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean deallocated arrays for ',kt,maskno
+          IF(lwp .AND. ln_diaregmean_ascii) CALL FLUSH(numdct_reg_txt)
+          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean flushed region mean text for ',kt,maskno
+      END DO
+
+      IF(lwp .AND. verbose) THEN                   ! Control print
+         WRITE(numout,*) 'dia_regmean ready to start iom_put'
+         CALL FLUSH(numout)
+      ENDIF
+      
+      !With current field_def.xml and iodef.xml, these fields must be output, so set to dummy values if not required.
+      
+      IF ( ln_diaregmean_nc ) THEN
+      
+          zrmet_out(:,:,:) = 0
+          zrmet_val = 0
+          tmp_name_iom = ''
+
+          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean ready to start iom_put: ',trim(tmp_name)
+          
+      
+          DO jm = 1,n_regions_output
+            zrmet_val = zrmet_ave(jm)
+!            if (zrmet_val .LT. -1e16) zrmet_val = -1e16
+!            if (zrmet_val .GT. 1e16) zrmet_val = 1e16
+            if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
+            zrmet_out(:,:,jm) = zrmet_val
+          END DO      
+          tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_ave'))
+          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom)
+          CALL iom_put(trim(tmp_name_iom), zrmet_out )
+          zrmet_out(:,:,:) = 0
+          zrmet_val = 0
+          tmp_name_iom = ''
+      
+          DO jm = 1,n_regions_output
+            zrmet_val = zrmet_tot(jm)
+!            if (zrmet_val .LT. -1e16) zrmet_val = -1e16
+!            if (zrmet_val .GT. 1e16) zrmet_val = 1e16
+            if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
+            zrmet_out(:,:,jm) = zrmet_val
+          END DO
+          tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_tot'))
+          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom)          
+          CALL iom_put( trim(tmp_name_iom), zrmet_out )
+          zrmet_out(:,:,:) = 0
+          zrmet_val = 0
+          tmp_name_iom = ''
+      
+          DO jm = 1,n_regions_output
+            zrmet_val = zrmet_var(jm)
+!            if (zrmet_val .LT. -1e16) zrmet_val = -1e16
+!            if (zrmet_val .GT. 1e16) zrmet_val = 1e16
+            if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
+            zrmet_out(:,:,jm) = zrmet_val
+          END DO
+          tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_var'))
+          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom)          
+          CALL iom_put( trim(tmp_name_iom), zrmet_out )
+          zrmet_out(:,:,:) = 0
+          zrmet_val = 0
+          tmp_name_iom = ''
+      
+          DO jm = 1,n_regions_output
+            zrmet_val = zrmet_cnt(jm)
+!            if (zrmet_val .LT. -1e16) zrmet_val = -1e16
+!            if (zrmet_val .GT. 1e16) zrmet_val = 1e16
+            if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
+            zrmet_out(:,:,jm) = zrmet_val
+          END DO
+          tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_cnt'))
+          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom)          
+          CALL iom_put( trim(tmp_name_iom), zrmet_out )
+          zrmet_out(:,:,:) = 0
+          zrmet_val = 0
+          tmp_name_iom = ''
+      
+          DO jm = 1,n_regions_output
+            zrmet_val = zrmet_reg_id(jm)
+!            if (zrmet_val .LT. -1e16) zrmet_val = -1e16
+!            if (zrmet_val .GT. 1e16) zrmet_val = 1e16
+            if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
+            zrmet_out(:,:,jm) = zrmet_val
+          END DO
+          tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_reg_id'))
+          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom)          
+          CALL iom_put( trim(tmp_name_iom), zrmet_out )
+          zrmet_out(:,:,:) = 0
+          zrmet_val = 0
+          tmp_name_iom = ''
+      
+          DO jm = 1,n_regions_output
+            zrmet_val = zrmet_mask_id(jm)
+!            if (zrmet_val .LT. -1e16) zrmet_val = -1e16
+!            if (zrmet_val .GT. 1e16) zrmet_val = 1e16
+            if (zrmet_val .NE. zrmet_val) zrmet_val = 1e20
+            zrmet_out(:,:,jm) = zrmet_val
+          END DO
+          tmp_name_iom =  trim(trim("reg_") // trim(tmp_name) // trim('_mask_id'))
+          IF(lwp .AND. verbose) WRITE(numout,*) 'dia_regmean iom_put tmp_name_iom : ',trim(tmp_name_iom)          
+          CALL iom_put( trim(tmp_name_iom), zrmet_out )
+          zrmet_out(:,:,:) = 0
+          zrmet_val = 0
+          tmp_name_iom = ''
+      ELSE
+        
+          ALLOCATE( dummy_zrmet(jpi,jpj,n_regions_output),  STAT= ierr )
+            IF( ierr /= 0 )   CALL ctl_stop( 'dia_wri_region_mean: failed to allocate dummy_zrmet array' )
+
+          DO jm = 1,n_regions_output
+              dummy_zrmet(:,:,jm) =     real(jm,wp)
+          END DO
+
+          DO jm = 1,9
+              CALL iom_put( trim(trim("reg_") // trim(tmp_name) // trim('_ave')), dummy_zrmet )
+              CALL iom_put( trim(trim("reg_") // trim(tmp_name) // trim('_tot')), dummy_zrmet )
+              CALL iom_put( trim(trim("reg_") // trim(tmp_name) // trim('_var')), dummy_zrmet )
+              CALL iom_put( trim(trim("reg_") // trim(tmp_name) // trim('_cnt')), dummy_zrmet )
+              CALL iom_put( trim(trim("reg_") // trim(tmp_name) // trim('_reg_id')), dummy_zrmet )
+              CALL iom_put( trim(trim("reg_") // trim(tmp_name) // trim('_mask_id')), dummy_zrmet )
+          END DO
+    
+          DEALLOCATE( dummy_zrmet)
+      ENDIF
+      
+      DEALLOCATE(zrmet_ave,zrmet_tot,zrmet_var,zrmet_cnt,zrmet_mask_id,zrmet_reg_id,zrmet_min,zrmet_max,zrmet_out)
+
+      IF(lwp .AND. verbose) THEN                   ! Control print
+         WRITE(numout,*) 
+         WRITE(numout,*) 'dia_wri_region_mean finished for ', trim(tmp_name)
+         WRITE(numout,*) 
+         CALL FLUSH(numout)
+      ENDIF
+      
+   END SUBROUTINE dia_wri_region_mean
+
+
+   !!======================================================================
+END MODULE diaregmean
Index: branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90	(revision 11635)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/DIA/diatmb.F90	(revision 11639)
@@ -11,4 +11,5 @@
    USE iom             ! I/0 library
    USE wrk_nemo        ! working arrays
+   USE diaregmean
 #if defined key_fabm
    USE trc, ONLY: trn
Index: branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 11635)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90	(revision 11639)
@@ -48,4 +48,6 @@
    USE dia25h          ! 25h Mean output
    USE diaopfoam       ! Diaopfoam output
+   USE diaregmean      ! regionalmean
+   USE diapea          ! pea
    USE iom
    USE ioipsl
@@ -406,4 +408,10 @@
       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
       !
Index: branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 11635)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 11639)
@@ -56,4 +56,6 @@
    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,4 +108,87 @@
       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,4 +312,11 @@
       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
Index: branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/OBS/obs_level_search.h90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/OBS/obs_level_search.h90	(revision 11635)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/OBS/obs_level_search.h90	(revision 11639)
@@ -21,4 +21,5 @@
       !!        !  2006-10  (A. Weaver) Cleanup
       !!        !  2008-10  (K. Mogensen) Remove assumptions on grid.
+      !!        !  2019-07  (R. Renshaw) Don't exceed model bottom.
       !!----------------------------------------------------------------------
 
@@ -46,5 +47,5 @@
             IF ( pgrddep(jk) >= pobsdep(ji) ) EXIT depk
          END DO depk
-         kobsk(ji) = jk
+         kobsk(ji) = MIN( jk, kgrd )
       END DO
 
Index: branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_z1d.h90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_z1d.h90	(revision 11635)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/OBS/obsinter_z1d.h90	(revision 11639)
@@ -26,4 +26,5 @@
       !!      ! 06-10 (A. Weaver) Cleanup
       !!      ! 07-01 (K. Mogensen) Use profile rather than single level
+      !!      ! 19-07 (R. Renshaw) Avoid division by zero
       !!---------------------------------------------------------------------
 
@@ -62,7 +63,12 @@
          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,4 +94,9 @@
 
          ENDIF
+
+         ELSE  ! take value directly from the higher model level
+            pobs(jdep)  = pobsk(kkco(jdep)-1)
+         ENDIF
+
       END DO
 
Index: branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 11635)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/nemogcm.F90	(revision 11639)
@@ -88,4 +88,6 @@
    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,4 +508,6 @@
       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,4 +647,5 @@
       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
 
       !
Index: branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/TOP_SRC/trcbdy.F90
===================================================================
--- branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/TOP_SRC/trcbdy.F90	(revision 11635)
+++ branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/TOP_SRC/trcbdy.F90	(revision 11639)
@@ -109,9 +109,14 @@
             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 )
