source: branches/UKMO/CO6_shelfclimate_fabm_noos/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 11264

Last change on this file since 11264 was 11264, checked in by hadjt, 15 months ago

Updated the maximum size of the NOOS cross-section size (from 150-275) in DIA/diadct.F90

  • Property svn:executable set to *
File size: 133.4 KB
Line 
1MODULE diadct
2
3  !!=====================================================================
4  !!                       ***  MODULE  diadct  ***
5  !! Ocean diagnostics: Compute the transport trough a sec.
6  !!===============================================================
7  !! History :
8  !!
9  !!         original  : 02/99 (Y Drillet)
10  !!         addition  : 10/01 (Y Drillet, R Bourdalle Badie)
11  !!                   : 10/05 (M Laborie) F90
12  !!         addition  : 04/07 (G Garric) Ice sections
13  !!         bugfix    : 04/07 (C Bricaud) test on sec%nb_point
14  !!                                      initialisation of ztransp1,ztransp2,...
15  !!         nemo_v_3_4: 09/2011 (C Bricaud)
16  !!
17  !!
18  !!----------------------------------------------------------------------
19#if defined key_diadct
20  !!----------------------------------------------------------------------
21  !!   'key_diadct' :
22  !!----------------------------------------------------------------------
23  !!----------------------------------------------------------------------
24  !!   dia_dct      :  Compute the transport through a sec.
25  !!   dia_dct_init :  Read namelist.
26  !!   readsec      :  Read sections description and pathway
27  !!   removepoints :  Remove points which are common to 2 procs
28  !!   transport    :  Compute transport for each sections
29  !!   dia_dct_wri  :  Write tranports results in ascii files
30  !!   interp       :  Compute temperature/salinity/density at U-point or V-point
31  !!   
32  !!----------------------------------------------------------------------
33  !! * Modules used
34  USE oce             ! ocean dynamics and tracers
35  USE dom_oce         ! ocean space and time domain
36  USE phycst          ! physical constants
37  USE in_out_manager  ! I/O units
38  USE iom             ! I/0 library
39  USE daymod          ! calendar
40  USE dianam          ! build name of file
41  USE lib_mpp         ! distributed memory computing library
42#if defined key_lim2
43  USE ice_2
44#endif
45#if defined key_lim3
46  USE ice
47#endif
48  USE domvvl
49  USE timing          ! preformance summary
50  USE wrk_nemo        ! working arrays
51 
52
53  IMPLICIT NONE
54  PRIVATE
55
56  !! * Routine accessibility
57  PUBLIC   dia_dct      ! routine called by step.F90
58  PUBLIC   dia_dct_init ! routine called by opa.F90
59  PUBLIC   diadct_alloc ! routine called by nemo_init in nemogcm.F90
60  PRIVATE  readsec
61  PRIVATE  removepoints
62  PRIVATE  transport
63  PRIVATE  dia_dct_wri
64  PRIVATE  dia_dct_wri_NOOS
65
66#include "domzgr_substitute.h90"
67
68  !! * Shared module variables
69  LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .TRUE.   !: model-data diagnostics flag
70  LOGICAL, PUBLIC ::   ln_dct_calc_noos_25h !: Calcuate noos 25 h means
71  LOGICAL, PUBLIC ::   ln_dct_calc_noos_hr !: Calcuate noos hourly means
72  ! JT
73  LOGICAL, PUBLIC ::   ln_dct_iom_cont !: Use IOM Output?
74  LOGICAL, PUBLIC ::   ln_dct_ascii    !: Output ascii or binary
75  LOGICAL, PUBLIC ::   ln_dct_h        !: Output hourly instantaneous or mean values
76  ! JT
77
78  !! * Module variables
79  INTEGER :: nn_dct        ! Frequency of computation
80  INTEGER :: nn_dctwri     ! Frequency of output
81  INTEGER :: nn_secdebug   ! Number of the section to debug
82!  INTEGER :: nn_dct_h    = 1     ! Frequency of computation for NOOS hourly files
83  INTEGER :: nn_dct_h      ! Frequency of computation for NOOS hourly files
84  INTEGER :: nn_dctwri_h   ! Frequency of output for NOOS hourly files
85   
86  INTEGER, PARAMETER :: nb_class_max  = 11   ! maximum number of classes, i.e. depth levels or density classes
87  ! JT INTEGER, PARAMETER :: nb_sec_max    = 30   ! maximum number of sections
88  INTEGER, PARAMETER :: nb_sec_max    = 80   !50 maximum number of sections
89  !JT INTEGER, PARAMETER :: nb_point_max  = 150  !375 maximum number of points in a single section
90  INTEGER, PARAMETER :: nb_point_max  = 275  !375 maximum number of points in a single section
91
92
93
94  INTEGER, PARAMETER :: nb_type       = 14   ! types of calculations, i.e. pos transport, neg transport, heat transport, salt transport
95  INTEGER, PARAMETER :: nb_3d_vars    = 5
96  INTEGER, PARAMETER :: nb_2d_vars    = 2
97  INTEGER            :: nb_sec 
98
99  TYPE POINT_SECTION
100     INTEGER :: I,J
101  END TYPE POINT_SECTION
102
103  TYPE COORD_SECTION
104     REAL(wp) :: lon,lat
105  END TYPE COORD_SECTION
106
107  TYPE SECTION
108     CHARACTER(len=60)                            :: name              ! name of the sec
109     LOGICAL                                      :: llstrpond         ! true if you want the computation of salt and
110                                                                       ! heat transports
111     LOGICAL                                      :: ll_ice_section    ! ice surface and ice volume computation
112     LOGICAL                                      :: ll_date_line      ! = T if the section crosses the date-line
113     TYPE(COORD_SECTION), DIMENSION(2)            :: coordSec          ! longitude and latitude of the extremities of the sec
114     INTEGER                                      :: nb_class          ! number of boundaries for density classes
115     INTEGER, DIMENSION(nb_point_max)             :: direction         ! vector direction of the point in the section
116     CHARACTER(len=40),DIMENSION(nb_class_max)    :: classname         ! characteristics of the class
117     REAL(wp), DIMENSION(nb_class_max)            :: zsigi           ,&! in-situ   density classes    (99 if you don't want)
118                                                     zsigp           ,&! potential density classes    (99 if you don't want)
119                                                     zsal            ,&! salinity classes   (99 if you don't want)
120                                                     ztem            ,&! temperature classes(99 if you don't want)
121                                                     zlay              ! level classes      (99 if you don't want)
122     REAL(wp), DIMENSION(nb_type,nb_class_max)        :: transport     ! transport output
123     REAL(wp), DIMENSION(nb_type,nb_class_max)        :: transport_h   ! transport output
124     REAL(wp)                                         :: slopeSection  ! slope of the section
125     INTEGER                                          :: nb_point      ! number of points in the section
126     TYPE(POINT_SECTION),DIMENSION(nb_point_max)      :: listPoint     ! list of points in the sections
127  END TYPE SECTION
128
129  TYPE(SECTION),DIMENSION(nb_sec_max) :: secs ! Array of sections
130
131  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  transports_3d
132  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d
133  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::  transports_3d_h
134  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  transports_2d_h
135  REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::  z_hr_output
136
137   !! $Id$
138CONTAINS
139
140  INTEGER FUNCTION diadct_alloc()
141     !!--------------------------------------------------------------------nb_sec--
142     !!                   ***  FUNCTION diadct_alloc  ***
143     !!----------------------------------------------------------------------
144     INTEGER :: ierr(2)
145     !!----------------------------------------------------------------------
146     !
147     
148     !JT not sure why this is in nemogcm.F90(?) rather than diadct_init...
149     !JT it would be good if the nb_sec_max and nb_point_max were controlled by name list variable.
150     
151     
152     ALLOCATE(transports_3d(nb_3d_vars,nb_sec_max,nb_point_max,jpk),   STAT=ierr(1) )
153     ALLOCATE(transports_2d(nb_2d_vars,nb_sec_max,nb_point_max)    ,   STAT=ierr(2) )
154     ALLOCATE(transports_3d_h(nb_3d_vars,nb_sec_max,nb_point_max,jpk), STAT=ierr(3) )
155     ALLOCATE(transports_2d_h(nb_2d_vars,nb_sec_max,nb_point_max)    , STAT=ierr(4) )
156     !JT ALLOCATE(z_hr_output(nb_sec_max,24,nb_class_max)                , STAT=ierr(5) )
157     ALLOCATE(z_hr_output(nb_sec_max,3,nb_class_max)                , STAT=ierr(5) )
158   
159     diadct_alloc = MAXVAL( ierr )
160     IF( diadct_alloc /= 0 )   CALL ctl_warn('diadct_alloc: failed to allocate arrays')
161     !
162  END FUNCTION diadct_alloc
163
164
165  SUBROUTINE dia_dct_init
166     !!---------------------------------------------------------------------
167     !!               ***  ROUTINE diadct  *** 
168     !!
169     !!  ** Purpose: Read the namelist parameters
170     !!              Open output files
171     !!
172     !!---------------------------------------------------------------------
173     !JT 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,nb_sec_max,nb_point_max
174     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
175     INTEGER           ::   ios,jsec        ! Local integer output status for namelist read
176     CHARACTER(len=3)  ::   jsec_str        ! name of the jsec
177
178     IF( nn_timing == 1 )   CALL timing_start('dia_dct_init')
179
180     REWIND( numnam_ref )              ! Namelist namdct in reference namelist : Diagnostic: transport through sections
181     READ  ( numnam_ref, namdct, IOSTAT = ios, ERR = 901)
182901  IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdct in reference namelist', lwp )
183
184     REWIND( numnam_cfg )              ! Namelist namdct in configuration namelist : Diagnostic: transport through sections
185     READ  ( numnam_cfg, namdct, IOSTAT = ios, ERR = 902 )
186902  IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdct in configuration namelist', lwp )
187     IF(lwm) WRITE ( numond, namdct )
188     
189
190     IF( ln_NOOS ) THEN
191
192        !Do calculation for daily, 25hourly mean every hour
193        nn_dct=3600./rdt         ! hard coded for NOOS transects, to give 25 hour means from hourly instantaneous values
194
195        !write out daily, 25hourly mean every day
196        nn_dctwri=86400./rdt
197
198
199        ! JT
200        !
201        !
202        !nn_dct_h=1       ! hard coded for NOOS transects, to give hourly data   
203        ! If you want hourly instantaneous values, you only do the calculation every 12 timesteps (if rdt = 300)
204        ! and output it every 12 time steps. For this, you set the ln_dct_h to be True, and it calcuates it automatically
205        ! if you want hourly mean values, set ln_dct_h to be False, and it will do the calculate every time step.
206        !
207        !SELECT CASE( ln_dct_h )
208        !   CASE(.TRUE.) 
209        !      nn_dct_h=3600./rdt
210        !   CASE(.FALSE.)
211        !      nn_dct_h=1
212        !END SELECT
213       
214        IF ( ln_dct_h ) THEN
215            nn_dct_h=3600./rdt
216        ELSE
217            nn_dct_h=1.
218        ENDIF   
219       
220        !JT write out hourly calculation every hour
221        nn_dctwri_h=3600./rdt
222     ENDIF
223
224     IF( lwp ) THEN
225        WRITE(numout,*) " "
226        WRITE(numout,*) "diadct_init: compute transports through sections "
227        WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
228        IF( ln_NOOS ) THEN
229           WRITE(numout,*) "       Calculate NOOS hourly output: ln_dct_calc_noos_hr = ",ln_dct_calc_noos_hr
230           WRITE(numout,*) "       Calculate NOOS 25 hour mean output: ln_dct_calc_noos_hr = ",ln_dct_calc_noos_25h
231           WRITE(numout,*) "       Use IOM Output: ln_dct_iom_cont = ",ln_dct_iom_cont
232           WRITE(numout,*) "       Output in ASCII (True) or Binary (False): ln_dct_ascii = ",ln_dct_ascii
233           WRITE(numout,*) "       Frequency of hourly computation - instantaneous (TRUE) or hourly mean (FALSE): ln_dct_h  = ",ln_dct_h
234
235           WRITE(numout,*) "       Frequency of computation hard coded to be every hour: nn_dct    = ",nn_dct
236           WRITE(numout,*) "       Frequency of write hard coded to average 25 instantaneous hour values: nn_dctwri = ",nn_dctwri
237           WRITE(numout,*) "       Frequency of hourly computation (timestep) : nn_dct_h  = ",nn_dct_h
238           WRITE(numout,*) "       Frequency of hourly computation Not hard coded to be every timestep, or : nn_dct_h  = ",nn_dct_h
239           WRITE(numout,*) "       Frequency of hourly write hard coded to every hour: nn_dctwri_h = ",nn_dctwri_h
240        ELSE
241           WRITE(numout,*) "       Frequency of computation: nn_dct    = ",nn_dct
242           WRITE(numout,*) "       Frequency of write:       nn_dctwri = ",nn_dctwri
243        ENDIF
244
245        IF      ( nn_secdebug .GE. 1 .AND. nn_secdebug .LE. nb_sec_max )THEN
246                                            WRITE(numout,*)"       Debug section number: ", nn_secdebug 
247        ELSE IF ( nn_secdebug ==  0 )THEN ; WRITE(numout,*)"       No section to debug"
248        ELSE IF ( nn_secdebug == -1 )THEN ; WRITE(numout,*)"       Debug all sections"
249        ELSE                              ; WRITE(numout,*)"       Wrong value for nn_secdebug : ",nn_secdebug
250        ENDIF
251
252        IF(nn_dct .GE. nn_dctwri .AND. MOD(nn_dct,nn_dctwri) .NE. 0)  &
253          &  CALL ctl_stop( 'diadct: nn_dct should be smaller and a multiple of nn_dctwri' )
254
255     ENDIF
256     
257     
258     IF ( ln_NOOS ) THEN
259        IF ( ln_dct_calc_noos_25h .or. ln_dct_calc_noos_hr ) CALL readsec
260     ENDIF
261
262     !open output file
263     IF( lwp ) THEN
264        IF( ln_NOOS ) THEN
265           WRITE(numout,*) "diadct_init: Open output files. ASCII? ",ln_dct_ascii
266           WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
267           IF  ( ln_dct_ascii ) THEN
268               if ( ln_dct_calc_noos_25h ) CALL ctl_opn( numdct_NOOS  ,'NOOS_transport'  , 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .TRUE. )
269               if ( ln_dct_calc_noos_hr )  CALL ctl_opn( numdct_NOOS_h,'NOOS_transport_h', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .TRUE. )
270           ELSE
271               if ( ln_dct_calc_noos_25h ) CALL ctl_opn( numdct_NOOS  ,'NOOS_transport_bin'  , 'REPLACE', 'UNFORMATTED', 'SEQUENTIAL', -1, numout,  .TRUE. )
272               if ( ln_dct_calc_noos_hr )  CALL ctl_opn( numdct_NOOS_h,'NOOS_transport_bin_h', 'REPLACE', 'UNFORMATTED', 'SEQUENTIAL', -1, numout,  .TRUE. )
273           ENDIF
274        ELSE
275           CALL ctl_opn( numdct_vol , 'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
276           CALL ctl_opn( numdct_temp, 'heat_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
277           CALL ctl_opn( numdct_sal , 'salt_transport'  , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout,  .FALSE. )
278        ENDIF
279     ENDIF
280
281     ! Initialise arrays to zero
282     transports_3d(:,:,:,:)   =0._wp
283     transports_2d(:,:,:)     =0._wp
284     transports_3d_h(:,:,:,:) =0._wp
285     transports_2d_h(:,:,:)   =0._wp
286     z_hr_output(:,:,:)       =0._wp
287
288     IF( nn_timing == 1 )   CALL timing_stop('dia_dct_init')
289
290     IF (ln_dct_iom_cont) THEN
291        IF( lwp ) THEN
292            WRITE(numout,*) " "
293            WRITE(numout,*) "diadct_init: using xios iom_put for output: field_def.xml and iodef.xml code"
294            WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
295            WRITE(numout,*) ""
296            WRITE(numout,*) "      field_def.xml"
297            WRITE(numout,*) "      ~~~~~~~~~~~~~"
298            WRITE(numout,*) ""
299            WRITE(numout,*) ""
300     
301            WRITE(numout,*)  '      <field_group id="noos_cross_section" domain_ref="1point" axis_ref="noos" operation="average">'
302
303            DO jsec=1,nb_sec
304               WRITE (jsec_str, "(I3.3)") jsec
305               
306               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" />'
307               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" />'
308               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" />'
309               
310            ENDDO
311           
312            WRITE(numout,*)  '      </field_group>'
313           
314            WRITE(numout,*) ""
315            WRITE(numout,*) ""
316            WRITE(numout,*) "      iodef.xml"
317            WRITE(numout,*) "      ~~~~~~~~~"
318            WRITE(numout,*) ""
319            WRITE(numout,*) ""
320           
321            WRITE(numout,*)  '      <file_group id="1d" output_freq="1d" output_level="10" enabled=".TRUE.">'
322            WRITE(numout,*) ""
323            WRITE(numout,*)  '          <file id="noos_cross_section" name="NOOS_transport">'
324            DO jsec=1,nb_sec
325               WRITE (jsec_str, "(I3.3)") jsec
326               
327               WRITE(numout,*)  '              <field field_ref="noos_'//jsec_str//'_trans" />'
328               WRITE(numout,*)  '              <field field_ref="noos_'//jsec_str//'_heat" />'
329               WRITE(numout,*)  '              <field field_ref="noos_'//jsec_str//'_salt" />'
330               
331            ENDDO
332            WRITE(numout,*)  '          </file>'
333            WRITE(numout,*) ""
334            WRITE(numout,*)  '      </file_group>'
335           
336            WRITE(numout,*) ""
337            WRITE(numout,*) ""
338            WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
339            WRITE(numout,*) ""
340       
341        ENDIF
342     ENDIF
343
344
345     !
346  END SUBROUTINE dia_dct_init
347 
348 
349  SUBROUTINE dia_dct(kt)
350     !!---------------------------------------------------------------------
351     !!               ***  ROUTINE diadct  *** 
352     !!
353     !!  Purpose :: Compute section transports and write it in numdct files
354     !! 
355     !!  Method  :: All arrays initialised to zero in dct_init
356     !!             Each nn_dct time step call subroutine 'transports' for
357     !!               each section to sum the transports.
358     !!             Each nn_dctwri time step:
359     !!               Divide the arrays by the number of summations to gain
360     !!               an average value
361     !!               Call dia_dct_sum to sum relevant grid boxes to obtain
362     !!               totals for each class (density, depth, temp or sal)
363     !!               Call dia_dct_wri to write the transports into file
364     !!               Reinitialise all relevant arrays to zero
365     !!---------------------------------------------------------------------
366     !! * Arguments
367     INTEGER,INTENT(IN)        ::kt
368
369     !! * Local variables
370     INTEGER             :: jsec,            &! loop on sections
371                            itotal            ! nb_sec_max*nb_type*nb_class_max
372     LOGICAL             :: lldebug =.FALSE.  ! debug a section 
373
374     
375     INTEGER , DIMENSION(1)             :: ish      ! tmp array for mpp_sum
376     INTEGER , DIMENSION(3)             :: ish2     !   "
377     REAL(wp), POINTER, DIMENSION(:)    :: zwork    !   "
378     REAL(wp), POINTER, DIMENSION(:,:,:):: zsum     !   "
379     !!---------------------------------------------------------------------   
380     
381     
382     
383     IF( nn_timing == 1 )   CALL timing_start('dia_dct')
384
385     IF( lk_mpp ) THEN
386        itotal = nb_sec_max*nb_type*nb_class_max
387        CALL wrk_alloc( itotal                          , zwork ) 
388        CALL wrk_alloc( nb_sec_max,nb_type,nb_class_max , zsum  )
389     ENDIF   
390 
391     ! Initialise arrays
392     zwork(:) = 0.0 
393     zsum(:,:,:) = 0.0
394
395     IF( lwp .AND. kt==nit000+nn_dct-1 ) THEN
396         WRITE(numout,*) " "
397         WRITE(numout,*) "diadct: compute transport"
398         WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~~~~~"
399         WRITE(numout,*) "nb_sec = ",nb_sec
400         WRITE(numout,*) "nn_dct = ",nn_dct
401         WRITE(numout,*) "ln_NOOS = ",ln_NOOS
402         WRITE(numout,*) "nb_sec = ",nb_sec
403         WRITE(numout,*) "nb_sec_max = ",nb_sec_max
404         WRITE(numout,*) "nb_type = ",nb_type
405         WRITE(numout,*) "nb_class_max = ",nb_class_max
406     ENDIF
407   
408   
409    IF ( ln_dct_calc_noos_25h ) THEN
410 
411        ! Compute transport and write only at nn_dctwri
412        IF ( MOD(kt,nn_dct)==0 .or. &               ! compute transport every nn_dct time steps
413            (ln_NOOS .and. kt==nn_it000 ) )  THEN   ! also include first time step when calculating NOOS 25 hour averages
414           
415           
416
417            DO jsec=1,nb_sec
418
419              lldebug=.FALSE.
420              IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND.  kt==nit000+nn_dct-1 .AND. lwp ) lldebug=.TRUE. 
421
422              !Compute transport through section 
423              CALL transport(secs(jsec),lldebug,jsec) 
424             
425
426            ENDDO
427               
428            IF( MOD(kt,nn_dctwri)==0 )THEN
429           
430           
431
432              IF( lwp .AND. kt==nit000+nn_dctwri-1 ) WRITE(numout,*)"      diadct: average and write at kt = ",kt
433
434
435              !JT   
436              !JT   
437              !JT   !! divide arrays by nn_dctwri/nn_dct to obtain average
438              !JT   transports_3d(:,:,:,:)=transports_3d(:,:,:,:)/(nn_dctwri/nn_dct)
439              !JT   transports_2d(:,:,:)  =transports_2d(:,:,:)  /(nn_dctwri/nn_dct)
440              !JT   
441              !JT   
442
443
444              !JT   
445              !JT Not 24 values, but 25! divide by ((nn_dctwri/nn_dct) +1)
446              !! divide arrays by nn_dctwri/nn_dct  to obtain average
447              transports_3d(:,:,:,:)= transports_3d(:,:,:,:)/((nn_dctwri/nn_dct)+1.)
448              transports_2d(:,:,:)  = transports_2d(:,:,:)  /((nn_dctwri/nn_dct)+1.)
449
450              ! Sum over each class
451              DO jsec=1,nb_sec
452                  CALL dia_dct_sum(secs(jsec),jsec)
453              ENDDO
454   
455              !Sum on all procs
456              IF( lk_mpp )THEN
457                  zsum(:,:,:)=0.0_wp
458                  ish(1)  =  nb_sec_max*nb_type*nb_class_max 
459                  ish2    = (/nb_sec_max,nb_type,nb_class_max/)
460                  DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport(:,:) ; ENDDO
461                  zwork(:)= RESHAPE(zsum(:,:,:), ish )
462                  CALL mpp_sum(zwork, ish(1))
463                  zsum(:,:,:)= RESHAPE(zwork,ish2)
464                  DO jsec=1,nb_sec ; secs(jsec)%transport(:,:) = zsum(jsec,:,:) ; ENDDO
465              ENDIF
466
467              !Write the transport
468              DO jsec=1,nb_sec
469
470                  IF( lwp .and. .not. ln_NOOS )CALL dia_dct_wri(kt,jsec,secs(jsec))
471                  !IF( lwp .and.       ln_NOOS )CALL dia_dct_wri_NOOS(kt,jsec,secs(jsec))   ! use NOOS specific formatting
472                  IF(  ln_NOOS )CALL dia_dct_wri_NOOS(kt,jsec,secs(jsec))   ! use NOOS specific formatting
473               
474               
475                  !nullify transports values after writing
476                  transports_3d(:,jsec,:,:)=0.0
477                  transports_2d(:,jsec,:  )=0.0
478                  secs(jsec)%transport(:,:)=0. 
479                 
480                 
481                  IF ( ln_NOOS ) CALL transport(secs(jsec),lldebug,jsec)  ! reinitialise for next 25 hour instantaneous average (overlapping values)
482
483
484
485              ENDDO
486
487            ENDIF
488
489        ENDIF
490       
491    ENDIF
492    IF ( ln_dct_calc_noos_hr ) THEN
493        IF ( MOD(kt,nn_dct_h)==0 ) THEN            ! compute transport every nn_dct_h time steps
494
495            DO jsec=1,nb_sec
496
497              lldebug=.FALSE.
498              IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND.  kt==nit000+nn_dct_h-1 .AND. lwp ) lldebug=.TRUE. 
499
500              !Compute transport through section 
501              CALL transport_h(secs(jsec),lldebug,jsec) 
502
503            ENDDO
504               
505            IF( MOD(kt,nn_dctwri_h)==0 )THEN
506
507              IF( lwp .AND. kt==nit000+nn_dctwri_h-1 )WRITE(numout,*)"      diadct: average and write hourly files at kt = ",kt         
508     
509              !! divide arrays by nn_dctwri/nn_dct to obtain average
510                !
511                ! JT - I think this is wrong. I think it is trying to sum over 25 hours, but only dividing by 24.
512                ! I think it might work for daily cycles, but not for monthly cycles,
513                !
514              transports_3d_h(:,:,:,:)=transports_3d_h(:,:,:,:)/(nn_dctwri_h/nn_dct_h)
515              transports_2d_h(:,:,:)  =transports_2d_h(:,:,:)  /(nn_dctwri_h/nn_dct_h)
516
517              ! Sum over each class
518              DO jsec=1,nb_sec
519                  CALL dia_dct_sum_h(secs(jsec),jsec)
520              ENDDO
521
522              !Sum on all procs
523              IF( lk_mpp )THEN
524                  ish(1)  =  nb_sec_max*nb_type*nb_class_max 
525                  ish2    = (/nb_sec_max,nb_type,nb_class_max/)
526                  DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport_h(:,:) ; ENDDO
527                  zwork(:)= RESHAPE(zsum(:,:,:), ish )
528                  CALL mpp_sum(zwork, ish(1))
529                  zsum(:,:,:)= RESHAPE(zwork,ish2)
530                  DO jsec=1,nb_sec ; secs(jsec)%transport_h(:,:) = zsum(jsec,:,:) ; ENDDO
531              ENDIF
532
533              !Write the transport
534              DO jsec=1,nb_sec
535
536                  IF( lwp .and. ln_NOOS ) THEN
537                    CALL dia_dct_wri_NOOS_h(kt/nn_dctwri_h,jsec,secs(jsec))   ! use NOOS specific formatting
538                  endif
539                  !nullify transports values after writing
540                  transports_3d_h(:,jsec,:,:)=0.0
541                  transports_2d_h(:,jsec,:)=0.0
542                  secs(jsec)%transport_h(:,:)=0.0
543                 
544                  ! for hourly mean or hourly instantaneous, you don't initialise! start with zero!
545                  !IF ( ln_NOOS ) CALL transport_h(secs(jsec),lldebug,jsec)  ! reinitialise for next 25 hour instantaneous average (overlapping values)
546
547              ENDDO
548
549            ENDIF
550
551        ENDIF   
552       
553     ENDIF
554
555     IF( lk_mpp )THEN
556        itotal = nb_sec_max*nb_type*nb_class_max
557        CALL wrk_dealloc( itotal                          , zwork ) 
558        CALL wrk_dealloc( nb_sec_max,nb_type,nb_class_max , zsum  )
559     ENDIF   
560
561     IF( nn_timing == 1 )   CALL timing_stop('dia_dct')
562     !
563  END SUBROUTINE dia_dct
564
565  SUBROUTINE readsec 
566     !!---------------------------------------------------------------------
567     !!               ***  ROUTINE readsec  ***
568     !!
569     !!  ** Purpose:
570     !!            Read a binary file(section_ijglobal.diadct)
571     !!            generated by the tools "NEMOGCM/TOOLS/SECTIONS_DIADCT"
572     !!
573     !!
574     !!---------------------------------------------------------------------
575     !! * Local variables
576     INTEGER :: iptglo , iptloc                               ! Global and local number of points for a section
577     INTEGER :: isec, iiglo, ijglo, iiloc, ijloc,iost,i1 ,i2  ! temporary  integer
578     INTEGER :: jsec, jpt                                     ! dummy loop indices
579                                                              ! heat/salt tranport is actived
580
581     INTEGER, DIMENSION(2) :: icoord 
582     CHARACTER(len=160)    :: clname                          !filename
583     CHARACTER(len=200)    :: cltmp
584     CHARACTER(len=200)    :: clformat                        !automatic format
585     TYPE(POINT_SECTION),DIMENSION(nb_point_max)  ::coordtemp !contains listpoints coordinates
586                                                              !read in the file
587     INTEGER, POINTER, DIMENSION(:) :: directemp              !contains listpoints directions
588                                                              !read in the files
589     LOGICAL :: llbon                                       ,&!local logical
590                lldebug                                       !debug the section
591     !!-------------------------------------------------------------------------------------
592     CALL wrk_alloc( nb_point_max, directemp )
593
594     !open input file
595     !---------------
596     !write(numout,*) 'dct low-level pre open: little endian '
597     !OPEN(UNIT=107,FILE='section_ijglobal.diadct', FORM='UNFORMATTED', ACCESS='SEQUENTIAL', STATUS='OLD',convert='LITTLE_ENDIAN')
598     
599     write(numout,*) 'dct low-level pre open: big endian :',nproc,narea
600     OPEN(UNIT=107,FILE='section_ijglobal.diadct', FORM='UNFORMATTED', ACCESS='SEQUENTIAL', STATUS='OLD',convert='BIG_ENDIAN')
601     
602     !write(numout,*) 'dct low-level pre open: SWAP '
603     !OPEN(UNIT=107,FILE='section_ijglobal.diadct', FORM='UNFORMATTED', ACCESS='SEQUENTIAL', STATUS='OLD',convert='SWAP')
604     
605     !write(numout,*) 'dct low-level pre open: NATIVE '
606     !OPEN(UNIT=107,FILE='section_ijglobal.diadct', FORM='UNFORMATTED', ACCESS='SEQUENTIAL', STATUS='OLD',convert='NATIVE')
607     
608     READ(107) isec
609     CLOSE(107)
610     
611     CALL ctl_opn( numdct_in, 'section_ijglobal.diadct', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, .TRUE. )
612 
613     
614     !---------------
615     !Read input file
616     !---------------
617     
618     DO jsec=1,nb_sec_max      !loop on the nb_sec sections
619
620        IF ( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) ) &
621           & WRITE(numout,*)'debugging for section number: ',jsec 
622
623        !initialization
624        !---------------
625        secs(jsec)%name=''
626        secs(jsec)%llstrpond      = .FALSE.
627        secs(jsec)%ll_ice_section = .FALSE.
628        secs(jsec)%ll_date_line   = .FALSE.
629        secs(jsec)%nb_class       = 0
630        secs(jsec)%zsigi          = 99._wp
631        secs(jsec)%zsigp          = 99._wp
632        secs(jsec)%zsal           = 99._wp
633        secs(jsec)%ztem           = 99._wp
634        secs(jsec)%zlay           = 99._wp
635        secs(jsec)%transport      =  0._wp
636        secs(jsec)%transport_h    =  0._wp
637        secs(jsec)%nb_point       = 0
638
639        !read section's number / name / computing choices / classes / slopeSection / points number
640        !-----------------------------------------------------------------------------------------
641       
642        READ(numdct_in,iostat=iost) isec
643        IF (iost .NE. 0 ) then
644          write(numout,*) 'unable to read section_ijglobal.diadct. iost = ',iost
645          EXIT !end of file
646        ENDIF
647       
648       
649        WRITE(cltmp,'(a,i4.4,a,i4.4)')'diadct: read sections : Problem of section number: isec= ',isec,' and jsec= ',jsec
650       
651       
652        IF( jsec .NE. isec )  CALL ctl_stop( cltmp )
653
654        READ(numdct_in)secs(jsec)%name
655        READ(numdct_in)secs(jsec)%llstrpond
656        READ(numdct_in)secs(jsec)%ll_ice_section
657        READ(numdct_in)secs(jsec)%ll_date_line
658        READ(numdct_in)secs(jsec)%coordSec
659        READ(numdct_in)secs(jsec)%nb_class
660        READ(numdct_in)secs(jsec)%zsigi
661        READ(numdct_in)secs(jsec)%zsigp
662        READ(numdct_in)secs(jsec)%zsal
663        READ(numdct_in)secs(jsec)%ztem
664        READ(numdct_in)secs(jsec)%zlay
665        READ(numdct_in)secs(jsec)%slopeSection
666        READ(numdct_in)iptglo
667
668        !debug
669        !-----
670
671        IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
672         
673            WRITE(clformat,'(a,i2,a)') '(A40,', nb_class_max,'(f8.3,1X))' 
674
675            WRITE(numout,*)       "   Section name :                       ",TRIM(secs(jsec)%name)
676            WRITE(numout,*)       "      Compute temperature and salinity transports ? ",secs(jsec)%llstrpond
677            WRITE(numout,*)       "      Compute ice transport ?           ",secs(jsec)%ll_ice_section
678            WRITE(numout,*)       "      Section crosses date-line ?       ",secs(jsec)%ll_date_line
679            WRITE(numout,*)       "      Slope section :                   ",secs(jsec)%slopeSection
680            WRITE(numout,*)       "      Number of points in the section:  ",iptglo
681            WRITE(numout,*)       "      Number of classes                 ",secs(jsec)%nb_class
682            WRITE(numout,clformat)"      Insitu density classes :          ",secs(jsec)%zsigi
683            WRITE(numout,clformat)"      Potential density classes :       ",secs(jsec)%zsigp
684            WRITE(numout,clformat)"      Salinity classes :                ",secs(jsec)%zsal
685            WRITE(numout,clformat)"      Temperature classes :             ",secs(jsec)%ztem
686            WRITE(numout,clformat)"      Depth classes :                   ",secs(jsec)%zlay
687        ENDIF     
688       
689
690        IF( iptglo .NE. 0 )THEN
691             
692           !read points'coordinates and directions
693           !--------------------------------------
694           coordtemp(:) = POINT_SECTION(0,0) !list of points read
695           directemp(:) = 0                  !value of directions of each points
696           DO jpt=1,iptglo
697              READ(numdct_in)i1,i2
698              coordtemp(jpt)%I = i1 
699              coordtemp(jpt)%J = i2
700           ENDDO
701           READ(numdct_in)directemp(1:iptglo)
702   
703           !debug
704           !-----
705           IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
706              WRITE(numout,*)"      List of points in global domain:"
707              DO jpt=1,iptglo
708                 WRITE(numout,*)'        # I J ',jpt,coordtemp(jpt),directemp(jpt)
709              ENDDO                 
710           ENDIF
711 
712           !Now each proc selects only points that are in its domain:
713           !--------------------------------------------------------
714           iptloc = 0                    !initialize number of points selected
715           DO jpt=1,iptglo               !loop on listpoint read in the file
716                   
717              iiglo=coordtemp(jpt)%I          ! global coordinates of the point
718              ijglo=coordtemp(jpt)%J          !  "
719
720              IF( iiglo==jpidta .AND. nimpp==1 ) iiglo = 2
721             
722              iiloc=iiglo-jpizoom+1-nimpp+1   ! local coordinates of the point
723              ijloc=ijglo-jpjzoom+1-njmpp+1   !  "
724
725              !verify if the point is on the local domain:(1,nlei)*(1,nlej)
726              IF( iiloc .GE. 1 .AND. iiloc .LE. nlei .AND. &
727              ijloc .GE. 1 .AND. ijloc .LE. nlej       )THEN
728
729                  WRITE(*,*)"diadct readsec: assigned proc!",narea,nproc,jpt
730                 
731                 
732                 iptloc = iptloc + 1                                                 ! count local points
733                 secs(jsec)%listPoint(iptloc) = POINT_SECTION(mi0(iiglo),mj0(ijglo)) ! store local coordinates
734                 secs(jsec)%direction(iptloc) = directemp(jpt)                       ! store local direction
735              ENDIF
736
737           ENDDO
738     
739           secs(jsec)%nb_point=iptloc !store number of section's points
740
741           !debug
742           !-----
743           IF(   lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
744              WRITE(numout,*)"      List of points selected by the proc:"
745              DO jpt = 1,iptloc
746                 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
747                 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1
748                 WRITE(numout,*)'         # I J : ',iiglo,ijglo
749              ENDDO
750           ENDIF
751
752           IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN
753              DO jpt = 1,iptloc
754                 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
755                 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1
756              ENDDO
757           ENDIF
758
759           !remove redundant points between processors
760           !------------------------------------------
761           lldebug = .FALSE. ; IF ( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND. lwp ) lldebug = .TRUE.
762           IF( iptloc .NE. 0 )THEN
763              CALL removepoints(secs(jsec),'I','top_list',lldebug)
764              CALL removepoints(secs(jsec),'I','bot_list',lldebug)
765              CALL removepoints(secs(jsec),'J','top_list',lldebug)
766              CALL removepoints(secs(jsec),'J','bot_list',lldebug)
767           ENDIF
768           IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN
769              DO jpt = 1,secs(jsec)%nb_point
770                 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
771                 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1
772              ENDDO
773           ENDIF
774
775           !debug
776           !-----
777           IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
778              WRITE(numout,*)"      List of points after removepoints:"
779              iptloc = secs(jsec)%nb_point
780              DO jpt = 1,iptloc
781                 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
782                 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1
783                 WRITE(numout,*)'         # I J : ',iiglo,ijglo
784                 CALL FLUSH(numout)
785              ENDDO
786           ENDIF
787
788        ELSE  ! iptglo = 0
789           IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )&
790              WRITE(numout,*)'   No points for this section.'
791        ENDIF
792
793     ENDDO !end of the loop on jsec
794     
795     nb_sec = jsec-1   !number of section read in the file
796
797     IF( lwp )  WRITE(numout,*)'diadct: read sections: Finished readsec.'
798
799     CALL wrk_dealloc( nb_point_max, directemp )
800     !
801  END SUBROUTINE readsec
802
803  SUBROUTINE removepoints(sec,cdind,cdextr,ld_debug)
804     !!---------------------------------------------------------------------------
805     !!             *** function removepoints
806     !!
807     !!   ** Purpose :: Remove points which are common to 2 procs
808     !!
809     !----------------------------------------------------------------------------
810     !! * arguments
811     TYPE(SECTION),INTENT(INOUT) :: sec
812     CHARACTER(len=1),INTENT(IN) :: cdind   ! = 'I'/'J'
813     CHARACTER(len=8),INTENT(IN) :: cdextr  ! = 'top_list'/'bot_list'
814     LOGICAL,INTENT(IN)          :: ld_debug                     
815
816     !! * Local variables
817     INTEGER :: iextr         ,& !extremity of listpoint that we verify
818                iind          ,& !coord     of listpoint that we verify
819                itest         ,& !indice value of the side of the domain
820                                 !where points could be redundant
821                isgn          ,& ! isgn= 1 : scan listpoint from start to end
822                                 ! isgn=-1 : scan listpoint from end to start
823                istart,iend      !first and last points selected in listpoint
824     INTEGER :: jpoint           !loop on list points
825     INTEGER, POINTER, DIMENSION(:)   :: idirec !contains temporary sec%direction
826     INTEGER, POINTER, DIMENSION(:,:) :: icoord !contains temporary sec%listpoint
827     !----------------------------------------------------------------------------
828     CALL wrk_alloc(    nb_point_max, idirec )
829     CALL wrk_alloc( 2, nb_point_max, icoord )
830
831     IF( ld_debug )WRITE(numout,*)'      -------------------------'
832     IF( ld_debug )WRITE(numout,*)'      removepoints in listpoint'
833
834     !iextr=extremity of list_point that we verify
835     IF      ( cdextr=='bot_list' )THEN ; iextr=1            ; isgn=1
836     ELSE IF ( cdextr=='top_list' )THEN ; iextr=sec%nb_point ; isgn=-1
837     ELSE    ; CALL ctl_stop("removepoints :Wrong value for cdextr")
838     ENDIF
839 
840     !which coordinate shall we verify ?
841     IF      ( cdind=='I' )THEN   ; itest=nlei ; iind=1
842     ELSE IF ( cdind=='J' )THEN   ; itest=nlej ; iind=2
843     ELSE    ; CALL ctl_stop("removepoints :Wrong value for cdind") 
844     ENDIF
845
846     IF( ld_debug )THEN
847        WRITE(numout,*)'      case: coord/list extr/domain side'
848        WRITE(numout,*)'      ', cdind,' ',cdextr,' ',itest
849        WRITE(numout,*)'      Actual number of points: ',sec%nb_point
850     ENDIF
851
852     icoord(1,1:nb_point_max) = sec%listPoint%I
853     icoord(2,1:nb_point_max) = sec%listPoint%J
854     idirec                   = sec%direction
855     sec%listPoint            = POINT_SECTION(0,0)
856     sec%direction            = 0
857
858     jpoint=iextr+isgn
859     DO WHILE( jpoint .GE. 1 .AND. jpoint .LE. sec%nb_point )
860         IF( icoord( iind,jpoint-isgn ) == itest .AND. icoord( iind,jpoint ) == itest )THEN ; jpoint=jpoint+isgn
861         ELSE                                                                               ; EXIT
862         ENDIF
863     ENDDO 
864
865     IF( cdextr=='bot_list')THEN ; istart=jpoint-1 ; iend=sec%nb_point
866     ELSE                        ; istart=1        ; iend=jpoint+1
867     ENDIF
868
869     sec%listPoint(1:1+iend-istart)%I = icoord(1,istart:iend)
870     sec%listPoint(1:1+iend-istart)%J = icoord(2,istart:iend)
871     sec%direction(1:1+iend-istart)   = idirec(istart:iend)
872     sec%nb_point                     = iend-istart+1
873     
874     IF( ld_debug )THEN
875        WRITE(numout,*)'      Number of points after removepoints :',sec%nb_point
876        WRITE(numout,*)'      sec%direction after removepoints :',sec%direction(1:sec%nb_point)
877     ENDIF
878
879     CALL wrk_dealloc(    nb_point_max, idirec )
880     CALL wrk_dealloc( 2, nb_point_max, icoord )
881
882  END SUBROUTINE removepoints
883
884  SUBROUTINE transport(sec,ld_debug,jsec)
885     !!-------------------------------------------------------------------------------------------
886     !!                     ***  ROUTINE transport  ***
887     !!
888     !!  Purpose ::  Compute the transport for each point in a section
889     !!
890     !!  Method  ::  Loop over each segment, and each vertical level and add the transport
891     !!              Be aware :         
892     !!              One section is a sum of segments
893     !!              One segment is defined by 2 consecutive points in sec%listPoint
894     !!              All points of sec%listPoint are positioned on the F-point of the cell
895     !!
896     !!              There are two loops:                 
897     !!              loop on the segment between 2 nodes
898     !!              loop on the level jk
899     !!
900     !! ** Output: Arrays containing the volume,density,salinity,temperature etc
901     !!            transports for each point in a section, summed over each nn_dct.
902     !!
903     !!-------------------------------------------------------------------------------------------
904     !! * Arguments
905     TYPE(SECTION),INTENT(INOUT) :: sec
906     LOGICAL      ,INTENT(IN)    :: ld_debug
907     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section
908   
909     !! * Local variables
910     INTEGER             :: jk,jseg,jclass,    &!loop on level/segment/classes
911                            isgnu  , isgnv      !
912     REAL(wp):: zumid        , zvmid        ,  &!U/V velocity on a cell segment
913                zumid_ice    , zvmid_ice    ,  &!U/V ice velocity
914                zTnorm                          !transport of velocity through one cell's sides
915     REAL(wp):: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point
916
917     TYPE(POINT_SECTION) :: k
918     !!--------------------------------------------------------
919
920     IF( ld_debug )WRITE(numout,*)'      Compute transport (jsec,sec%nb_point,sec%slopeSection) : ', jsec,sec%nb_point,sec%slopeSection
921
922     !---------------------------!
923     !  COMPUTE TRANSPORT        !
924     !---------------------------!
925     IF(sec%nb_point .NE. 0)THEN
926
927        !----------------------------------------------------------------------------------------------------
928        !----------------------------------------------------------------------------------------------------
929        !----------------------------------------------------------------------------------------------------
930        !
931        !
932        ! ! ! ! JT 1/09/2018 - changing convention. Always direction + is toward left hand of section
933        !
934        !    Making sign of the velocities used to calculate the volume transport a function of direction, not slopesection
935        !    (isgnu, isgnv)
936        !   
937        !    They vary for each segment of the section.
938        !
939        !----------------------------------------------------------------------------------------------------
940        !----------------------------------------------------------------------------------------------------
941        !----------------------------------------------------------------------------------------------------
942        !Compute sign for velocities:
943        !
944        !convention:
945        !   non horizontal section: direction + is toward left hand of section
946        !       horizontal section: direction + is toward north of section
947        !
948        !
949        !       slopeSection < 0     slopeSection > 0       slopeSection=inf            slopeSection=0
950        !       ----------------      -----------------     ---------------             --------------
951        !
952        !   isgnv=1         direction +     
953        !  ______         _____             ______                                                   
954        !        |           //|            |                  |                         direction +   
955        !        | isgnu=1  // |            |isgnu=1           |isgnu=1                     /|\
956        !        |_______  //         ______|    \\            | ---\                        |
957        !               |             | isgnv=-1  \\ |         | ---/ direction +       ____________
958        !               |             |          __\\|         |                   
959        !               |             |     direction +        |                      isgnv=1                                 
960        !                                                     
961        !----------------------------------------------------------------------------------------------------
962        ! JT      isgnu = 1
963        ! JT      IF( sec%slopeSection .GT. 0 ) THEN  ; isgnv = -1
964        ! JT      ELSE                                ; isgnv =  1
965        ! JT      ENDIF
966        ! JT      IF( sec%slopeSection .GE. 9999. )     isgnv =  1
967
968        IF( ld_debug )write(numout,*)"sec%slopeSection isgnu isgnv ",sec%slopeSection,isgnu,isgnv
969
970        !--------------------------------------!
971        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  !
972        !--------------------------------------!
973        DO jseg=1,MAX(sec%nb_point-1,0)
974           
975           
976           !JT: Compute sign for velocities:
977           
978           !isgnu =  1
979           !isgnv =  1
980           !
981           ! JT changing sign of u and v is dependent on the direction of the section.
982           !isgnu =  1
983           !isgnv =  1
984           !SELECT CASE( sec%direction(jseg) )
985           !CASE(0)  ;   isgnv = -1
986           !CASE(3)  ;   isgnu = -1
987           !END SELECT
988           
989           
990           SELECT CASE( sec%direction(jseg) )
991           CASE(0) 
992              isgnu =  1
993              isgnv = -1
994           CASE(1)
995              isgnu =  1
996              isgnv =  1
997           CASE(2) 
998              isgnu =  1
999              isgnv =  1
1000           CASE(3) 
1001              isgnu = -1
1002              isgnv =  1
1003           END SELECT
1004           
1005           !-------------------------------------------------------------------------------------------
1006           ! Select the appropriate coordinate for computing the velocity of the segment
1007           ! Corrected by JT 01/09/2018 (#)
1008           !
1009           !                      CASE(0)                                    Case (2)
1010           !                      -------                                    --------
1011           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)     
1012           !      F(i,j)---------#V(i,j)-------F(i+1,j)                                 |
1013           !                     -------->                                              |
1014           !                                                                        |   |
1015           !                                                                        |   |
1016           !                      Case (3)                                          | U(i,j)
1017           !                      --------                                          |   |
1018           !                                                                        V   |
1019           !  listPoint(jseg+1) F(i,j+1)                                                |
1020           !                        |                                                   |
1021           !                        |                                                   |
1022           !                        |                                 listPoint(jseg+1) F(i,j-1)
1023           !                   ^    |                                           
1024           !                   |    |                                           
1025           !                   | U(i,j+1)                                           
1026           !                   |    |                                       Case(1)     
1027           !                   |    |                                       ------     
1028           !                        |                                           
1029           !                        |                 listPoint(jseg+1)             listPoint(jseg)                           
1030           !                        |                 F(i-1,j)----------#V(i-1,j) ------#f(i,j)                           
1031           ! listPoint(jseg)     F(i,j)                                 <-------
1032           !
1033           !-------------------------------------------------------------------------------------------
1034
1035           SELECT CASE( sec%direction(jseg) )
1036           CASE(0)  ;   k = sec%listPoint(jseg)
1037           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
1038           CASE(2)  ;   k = sec%listPoint(jseg)
1039           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
1040           END SELECT
1041
1042           !---------------------------|
1043           !     LOOP ON THE LEVEL     |
1044           !---------------------------|
1045           !Sum of the transport on the vertical
1046           DO jk=1,mbathy(k%I,k%J)
1047
1048              ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
1049              SELECT CASE( sec%direction(jseg) )
1050              CASE(0,1)
1051                 ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) )
1052                 zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) )
1053                 zrhop = interp(k%I,k%J,jk,'V',rhop)
1054                 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0)
1055                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1)
1056              CASE(2,3)
1057                 ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) )
1058                 zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) )
1059                 zrhop = interp(k%I,k%J,jk,'U',rhop)
1060                 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0)
1061                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1) 
1062              END SELECT
1063
1064              zfsdep= fsdept(k%I,k%J,jk)
1065 
1066              !compute velocity with the correct direction
1067              SELECT CASE( sec%direction(jseg) )
1068              CASE(0,1) 
1069                 zumid=0.
1070                 zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk)
1071              CASE(2,3)
1072                 zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk)
1073                 zvmid=0.
1074              END SELECT
1075
1076              !zTnorm=transport through one cell;
1077              !velocity* cell's length * cell's thickness
1078              zTnorm=zumid*e2u(k%I,k%J)*  fse3u(k%I,k%J,jk)+     &
1079                     zvmid*e1v(k%I,k%J)*  fse3v(k%I,k%J,jk)
1080
1081#if ! defined key_vvl
1082              !add transport due to free surface
1083              IF( jk==1 )THEN
1084                 zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + &
1085                                   zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk)
1086              ENDIF
1087#endif
1088              !COMPUTE TRANSPORT
1089
1090              transports_3d(1,jsec,jseg,jk) = transports_3d(1,jsec,jseg,jk) + zTnorm
1091 
1092              IF ( sec%llstrpond ) THEN
1093                 transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk)  + zTnorm * zrhoi
1094                 transports_3d(3,jsec,jseg,jk) = transports_3d(3,jsec,jseg,jk)  + zTnorm * zrhop
1095                 transports_3d(4,jsec,jseg,jk) = transports_3d(4,jsec,jseg,jk)  + zTnorm * 4.e+3_wp * (ztn+273.15) * 1026._wp
1096                 transports_3d(5,jsec,jseg,jk) = transports_3d(5,jsec,jseg,jk)  + zTnorm * 0.001 * zsn * 1026._wp
1097              ENDIF
1098
1099           ENDDO !end of loop on the level
1100
1101#if defined key_lim2 || defined key_lim3
1102
1103           !ICE CASE   
1104           !------------
1105           IF( sec%ll_ice_section )THEN
1106              SELECT CASE (sec%direction(jseg))
1107              CASE(0)
1108                 zumid_ice = 0
1109                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
1110              CASE(1)
1111                 zumid_ice = 0
1112                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
1113              CASE(2)
1114                 zvmid_ice = 0
1115                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
1116              CASE(3)
1117                 zvmid_ice = 0
1118                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
1119              END SELECT
1120   
1121              zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J)
1122   
1123              transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)*   &
1124                                   (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &
1125                                  *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) +  &
1126                                    hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) &
1127                                   +zice_vol_pos
1128              transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)*   &
1129                                    (1.0 -  frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &
1130                                   +zice_surf_pos
1131   
1132           ENDIF !end of ice case
1133#endif
1134 
1135        ENDDO !end of loop on the segment
1136
1137     ENDIF   !end of sec%nb_point =0 case
1138     !
1139  END SUBROUTINE transport
1140
1141  SUBROUTINE transport_h(sec,ld_debug,jsec)
1142     !!-------------------------------------------------------------------------------------------
1143     !!                     ***  ROUTINE hourly_transport  ***
1144     !!
1145     !!              Exactly as routine transport but for data summed at
1146     !!              each time step and averaged each hour
1147     !!
1148     !!  Purpose ::  Compute the transport for each point in a section
1149     !!
1150     !!  Method  ::  Loop over each segment, and each vertical level and add the transport
1151     !!              Be aware :         
1152     !!              One section is a sum of segments
1153     !!              One segment is defined by 2 consecutive points in sec%listPoint
1154     !!              All points of sec%listPoint are positioned on the F-point of the cell
1155     !!
1156     !!              There are two loops:                 
1157     !!              loop on the segment between 2 nodes
1158     !!              loop on the level jk
1159     !!
1160     !! ** Output: Arrays containing the volume,density,salinity,temperature etc
1161     !!            transports for each point in a section, summed over each nn_dct.
1162     !!
1163     !!-------------------------------------------------------------------------------------------
1164     !! * Arguments
1165     TYPE(SECTION),INTENT(INOUT) :: sec
1166     LOGICAL      ,INTENT(IN)    :: ld_debug
1167     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section
1168   
1169     !! * Local variables
1170     INTEGER             :: jk,jseg,jclass,    &!loop on level/segment/classes
1171                            isgnu  , isgnv      !
1172     REAL(wp):: zumid        , zvmid        ,  &!U/V velocity on a cell segment
1173                zumid_ice    , zvmid_ice    ,  &!U/V ice velocity
1174                zTnorm                          !transport of velocity through one cell's sides
1175     REAL(wp):: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point
1176
1177     TYPE(POINT_SECTION) :: k
1178     !!--------------------------------------------------------
1179
1180     IF( ld_debug )WRITE(numout,*)'      Compute transport'
1181
1182     !---------------------------!
1183     !  COMPUTE TRANSPORT        !
1184     !---------------------------!
1185     IF(sec%nb_point .NE. 0)THEN
1186
1187        !----------------------------------------------------------------------------------------------------
1188        !----------------------------------------------------------------------------------------------------
1189        !----------------------------------------------------------------------------------------------------
1190        !
1191        !
1192        ! ! ! ! JT 1/09/2018 - changing convention. Always direction + is toward left hand of section
1193        !
1194        !    Making sign of the velocities used to calculate the volume transport a function of direction, not slopesection
1195        !    (isgnu, isgnv)
1196        !   
1197        !    They vary for each segment of the section.
1198        !
1199        !----------------------------------------------------------------------------------------------------
1200        !----------------------------------------------------------------------------------------------------
1201        !----------------------------------------------------------------------------------------------------
1202        !Compute sign for velocities:
1203        !
1204        !convention:
1205        !   non horizontal section: direction + is toward left hand of section
1206        !       horizontal section: direction + is toward north of section
1207        !
1208        !
1209        !       slopeSection < 0     slopeSection > 0       slopeSection=inf            slopeSection=0
1210        !       ----------------      -----------------     ---------------             --------------
1211        !
1212        !   isgnv=1         direction +     
1213        !  ______         _____             ______                                                   
1214        !        |           //|            |                  |                         direction +   
1215        !        | isgnu=1  // |            |isgnu=1           |isgnu=1                     /|\
1216        !        |_______  //         ______|    \\            | ---\                        |
1217        !               |             | isgnv=-1  \\ |         | ---/ direction +       ____________
1218        !               |             |          __\\|         |                   
1219        !               |             |     direction +        |                      isgnv=1                                 
1220        !                                                     
1221        !----------------------------------------------------------------------------------------------------
1222        ! JT      isgnu = 1
1223        ! JT      IF( sec%slopeSection .GT. 0 ) THEN  ; isgnv = -1
1224        ! JT      ELSE                                ; isgnv =  1
1225        ! JT      ENDIF
1226        ! JT      IF( sec%slopeSection .GE. 9999. )     isgnv =  1
1227
1228        IF( ld_debug )write(numout,*)"isgnu isgnv ",isgnu,isgnv
1229
1230        !--------------------------------------!
1231        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  !
1232        !--------------------------------------!
1233        DO jseg=1,MAX(sec%nb_point-1,0)
1234           
1235           
1236           !JT: Compute sign for velocities:
1237           
1238           !isgnu =  1
1239           !isgnv =  1
1240           !
1241           ! JT changing sign of u and v is dependent on the direction of the section.
1242           !isgnu =  1
1243           !isgnv =  1
1244           !SELECT CASE( sec%direction(jseg) )
1245           !CASE(0)  ;   isgnv = -1
1246           !CASE(3)  ;   isgnu = -1
1247           !END SELECT
1248           
1249           
1250           SELECT CASE( sec%direction(jseg) )
1251           CASE(0) 
1252              isgnu =  1
1253              isgnv = -1
1254           CASE(1)
1255              isgnu =  1
1256              isgnv =  1
1257           CASE(2) 
1258              isgnu =  1
1259              isgnv =  1
1260           CASE(3) 
1261              isgnu = -1
1262              isgnv =  1
1263           END SELECT
1264           
1265           !-------------------------------------------------------------------------------------------
1266           ! Select the appropriate coordinate for computing the velocity of the segment
1267           ! Corrected by JT 01/09/2018 (#)
1268           !
1269           !                      CASE(0)                                    Case (2)
1270           !                      -------                                    --------
1271           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)     
1272           !      F(i,j)---------#V(i,j)-------F(i+1,j)                                 |
1273           !                     -------->                                              |
1274           !                                                                        |   |
1275           !                                                                        |   |
1276           !                      Case (3)                                          | U(i,j)
1277           !                      --------                                          |   |
1278           !                                                                        V   |
1279           !  listPoint(jseg+1) F(i,j+1)                                                |
1280           !                        |                                                   |
1281           !                        |                                                   |
1282           !                        |                                 listPoint(jseg+1) F(i,j-1)
1283           !                   ^    |                                           
1284           !                   |    |                                           
1285           !                   | U(i,j+1)                                           
1286           !                   |    |                                       Case(1)     
1287           !                   |    |                                       ------     
1288           !                        |                                           
1289           !                        |                 listPoint(jseg+1)             listPoint(jseg)                           
1290           !                        |                 F(i-1,j)----------#V(i-1,j) ------#f(i,j)                           
1291           ! listPoint(jseg)     F(i,j)                                 <-------
1292           !
1293           !-------------------------------------------------------------------------------------------
1294
1295           SELECT CASE( sec%direction(jseg) )
1296           CASE(0)  ;   k = sec%listPoint(jseg)
1297           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
1298           CASE(2)  ;   k = sec%listPoint(jseg)
1299           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
1300           END SELECT
1301
1302           !---------------------------|
1303           !     LOOP ON THE LEVEL     |
1304           !---------------------------|
1305           !Sum of the transport on the vertical
1306           DO jk=1,mbathy(k%I,k%J)
1307
1308              ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
1309              SELECT CASE( sec%direction(jseg) )
1310              CASE(0,1)
1311                 ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) )
1312                 zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) )
1313                 zrhop = interp(k%I,k%J,jk,'V',rhop)
1314                 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0)
1315                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1)
1316              CASE(2,3)
1317                 ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) )
1318                 zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) )
1319                 zrhop = interp(k%I,k%J,jk,'U',rhop)
1320                 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0)
1321                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1) 
1322              END SELECT
1323
1324              zfsdep= fsdept(k%I,k%J,jk)
1325 
1326              !compute velocity with the correct direction
1327              SELECT CASE( sec%direction(jseg) )
1328              CASE(0,1) 
1329                 zumid=0.
1330                 zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk)
1331              CASE(2,3)
1332                 zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk)
1333                 zvmid=0.
1334              END SELECT
1335
1336              !zTnorm=transport through one cell;
1337              !velocity* cell's length * cell's thickness
1338              zTnorm=zumid*e2u(k%I,k%J)*  fse3u(k%I,k%J,jk)+     &
1339                     zvmid*e1v(k%I,k%J)*  fse3v(k%I,k%J,jk)
1340
1341#if ! defined key_vvl
1342              !add transport due to free surface
1343              IF( jk==1 )THEN
1344                 zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + &
1345                                   zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk)
1346              ENDIF
1347#endif
1348              !COMPUTE TRANSPORT
1349
1350              transports_3d_h(1,jsec,jseg,jk) = transports_3d_h(1,jsec,jseg,jk) + zTnorm
1351 
1352              IF ( sec%llstrpond ) THEN
1353                 transports_3d_h(2,jsec,jseg,jk) = transports_3d_h(2,jsec,jseg,jk)  + zTnorm * zrhoi
1354                 transports_3d_h(3,jsec,jseg,jk) = transports_3d_h(3,jsec,jseg,jk)  + zTnorm * zrhop
1355                 transports_3d_h(4,jsec,jseg,jk) = transports_3d_h(4,jsec,jseg,jk)  + zTnorm * 4.e+3_wp * (ztn+273.15) * 1026._wp
1356                 transports_3d_h(5,jsec,jseg,jk) = transports_3d_h(5,jsec,jseg,jk)  + zTnorm * 0.001 * zsn * 1026._wp
1357              ENDIF
1358
1359           ENDDO !end of loop on the level
1360
1361#if defined key_lim2 || defined key_lim3
1362
1363           !ICE CASE   
1364           !------------
1365           IF( sec%ll_ice_section )THEN
1366              SELECT CASE (sec%direction(jseg))
1367              CASE(0)
1368                 zumid_ice = 0
1369                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
1370              CASE(1)
1371                 zumid_ice = 0
1372                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
1373              CASE(2)
1374                 zvmid_ice = 0
1375                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
1376              CASE(3)
1377                 zvmid_ice = 0
1378                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
1379              END SELECT
1380   
1381              zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J)
1382   
1383              transports_2d_h(1,jsec,jseg) = transports_2d_h(1,jsec,jseg) + (zTnorm)*   &
1384                                   (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &
1385                                  *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) +  &
1386                                    hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) &
1387                                   +zice_vol_pos
1388              transports_2d_h(2,jsec,jseg) = transports_2d_h(2,jsec,jseg) + (zTnorm)*   &
1389                                    (1.0 -  frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &
1390                                   +zice_surf_pos
1391   
1392           ENDIF !end of ice case
1393#endif
1394 
1395        ENDDO !end of loop on the segment
1396
1397     ENDIF   !end of sec%nb_point =0 case
1398     !
1399  END SUBROUTINE transport_h
1400 
1401  SUBROUTINE dia_dct_sum(sec,jsec)
1402     !!-------------------------------------------------------------
1403     !! Purpose: Average the transport over nn_dctwri time steps
1404     !! and sum over the density/salinity/temperature/depth classes
1405     !!
1406     !! Method:
1407     !!           Sum over relevant grid cells to obtain values
1408     !!           for each
1409     !!              There are several loops:                 
1410     !!              loop on the segment between 2 nodes
1411     !!              loop on the level jk
1412     !!              loop on the density/temperature/salinity/level classes
1413     !!              test on the density/temperature/salinity/level
1414     !!
1415     !!  ** Method  :Transport through a given section is equal to the sum of transports
1416     !!              computed on each proc.
1417     !!              On each proc,transport is equal to the sum of transport computed through
1418     !!              segments linking each point of sec%listPoint  with the next one.   
1419     !!
1420     !!-------------------------------------------------------------
1421     !! * arguments
1422     TYPE(SECTION),INTENT(INOUT) :: sec
1423     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section
1424
1425     TYPE(POINT_SECTION) :: k
1426     INTEGER  :: jk,jseg,jclass                        !loop on level/segment/classes
1427     REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point
1428     !!-------------------------------------------------------------
1429
1430     !! Sum the relevant segments to obtain values for each class
1431     IF(sec%nb_point .NE. 0)THEN   
1432
1433        !--------------------------------------!
1434        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  !
1435        !--------------------------------------!
1436        DO jseg=1,MAX(sec%nb_point-1,0)
1437           
1438           !-------------------------------------------------------------------------------------------
1439           ! Select the appropriate coordinate for computing the velocity of the segment
1440           !
1441           !                      CASE(0)                                    Case (2)
1442           !                      -------                                    --------
1443           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)     
1444           !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               |
1445           !                                                                            |
1446           !                                                                            |
1447           !                                                                            |
1448           !                      Case (3)                                            U(i,j)
1449           !                      --------                                              |
1450           !                                                                            |
1451           !  listPoint(jseg+1) F(i,j+1)                                                |
1452           !                        |                                                   |
1453           !                        |                                                   |
1454           !                        |                                 listPoint(jseg+1) F(i,j-1)
1455           !                        |                                           
1456           !                        |                                           
1457           !                     U(i,j+1)                                           
1458           !                        |                                       Case(1)     
1459           !                        |                                       ------     
1460           !                        |                                           
1461           !                        |                 listPoint(jseg+1)             listPoint(jseg)                           
1462           !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                           
1463           ! listPoint(jseg)     F(i,j)
1464           !
1465           !-------------------------------------------------------------------------------------------
1466
1467           SELECT CASE( sec%direction(jseg) )
1468           CASE(0)  ;   k = sec%listPoint(jseg)
1469           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
1470           CASE(2)  ;   k = sec%listPoint(jseg)
1471           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
1472           END SELECT
1473
1474           !---------------------------|
1475           !     LOOP ON THE LEVEL     |
1476           !---------------------------|
1477           !Sum of the transport on the vertical
1478           DO jk=1,mbathy(k%I,k%J)
1479
1480              ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
1481              SELECT CASE( sec%direction(jseg) )
1482              CASE(0,1)
1483                 ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) )
1484                 zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) )
1485                 zrhop = interp(k%I,k%J,jk,'V',rhop)
1486                 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0)
1487                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1)
1488              CASE(2,3)
1489                 ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) )
1490                 zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) )
1491                 zrhop = interp(k%I,k%J,jk,'U',rhop)
1492                 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0)
1493                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1) 
1494              END SELECT
1495
1496              zfsdep= fsdept(k%I,k%J,jk) 
1497 
1498              !-------------------------------
1499              !  LOOP ON THE DENSITY CLASSES |
1500              !-------------------------------
1501              !The computation is made for each density/temperature/salinity/depth class
1502              DO jclass=1,MAX(1,sec%nb_class-1)
1503 
1504                 !----------------------------------------------!
1505                 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL!
1506                 !----------------------------------------------!
1507 
1508                 IF ( (                                                    &
1509                    ((( zrhop .GE. (sec%zsigp(jclass)+1000.  )) .AND.      &
1510                    (   zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR.     &
1511                    ( sec%zsigp(jclass) .EQ. 99.)) .AND.                   &
1512
1513                    ((( zrhoi .GE. (sec%zsigi(jclass) + 1000.  )) .AND.    &
1514                    (   zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR.     &
1515                    ( sec%zsigi(jclass) .EQ. 99.)) .AND.                   &
1516
1517                    ((( zsn .GT. sec%zsal(jclass)) .AND.                   &
1518                    (   zsn .LE. sec%zsal(jclass+1))) .OR.                 &
1519                    ( sec%zsal(jclass) .EQ. 99.)) .AND.                    &
1520
1521                    ((( ztn .GE. sec%ztem(jclass)) .AND.                   &
1522                    (   ztn .LE. sec%ztem(jclass+1))) .OR.                 &
1523                    ( sec%ztem(jclass) .EQ.99.)) .AND.                     &
1524
1525                    ((( zfsdep .GE. sec%zlay(jclass)) .AND.                &
1526                    (   zfsdep .LE. sec%zlay(jclass+1))) .OR.              &
1527                    ( sec%zlay(jclass) .EQ. 99. ))                         &
1528                                                                   ))   THEN
1529
1530                    !SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS
1531                    !----------------------------------------------------------------------------
1532                    IF (transports_3d(1,jsec,jseg,jk) .GE. 0.0) THEN
1533                       sec%transport(1,jclass) = sec%transport(1,jclass)+transports_3d(1,jsec,jseg,jk)
1534                    ELSE
1535                       sec%transport(2,jclass) = sec%transport(2,jclass)+transports_3d(1,jsec,jseg,jk)
1536                    ENDIF
1537                    IF( sec%llstrpond )THEN
1538
1539                       IF( transports_3d(1,jsec,jseg,jk) .NE. 0._wp ) THEN
1540
1541                          IF (transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) .GE. 0.0 ) THEN
1542                             sec%transport(3,jclass) = sec%transport(3,jclass)+transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk)
1543                          ELSE
1544                             sec%transport(4,jclass) = sec%transport(4,jclass)+transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk)
1545                          ENDIF
1546
1547                          IF ( transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) .GE. 0.0 ) THEN
1548                             sec%transport(5,jclass) = sec%transport(5,jclass)+transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk)
1549                          ELSE
1550                             sec%transport(6,jclass) = sec%transport(6,jclass)+transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk)
1551                          ENDIF
1552
1553                       ENDIF
1554
1555                       IF ( transports_3d(4,jsec,jseg,jk) .GE. 0.0 ) THEN
1556                          sec%transport(7,jclass) = sec%transport(7,jclass)+transports_3d(4,jsec,jseg,jk)
1557                       ELSE
1558                          sec%transport(8,jclass) = sec%transport(8,jclass)+transports_3d(4,jsec,jseg,jk)
1559                       ENDIF
1560
1561                       IF ( transports_3d(5,jsec,jseg,jk) .GE. 0.0 ) THEN
1562                          sec%transport( 9,jclass) = sec%transport( 9,jclass)+transports_3d(5,jsec,jseg,jk)
1563                       ELSE
1564                          sec%transport(10,jclass) = sec%transport(10,jclass)+transports_3d(5,jsec,jseg,jk)
1565                       ENDIF
1566
1567                    ELSE
1568                       sec%transport( 3,jclass) = 0._wp
1569                       sec%transport( 4,jclass) = 0._wp
1570                       sec%transport( 5,jclass) = 0._wp
1571                       sec%transport( 6,jclass) = 0._wp
1572                       sec%transport( 7,jclass) = 0._wp
1573                       sec%transport( 8,jclass) = 0._wp
1574                       sec%transport( 9,jclass) = 0._wp
1575                       sec%transport(10,jclass) = 0._wp
1576                    ENDIF
1577
1578                 ENDIF ! end of test if point is in class
1579   
1580              ENDDO ! end of loop on the classes
1581
1582           ENDDO ! loop over jk
1583
1584#if defined key_lim2 || defined key_lim3
1585
1586           !ICE CASE   
1587           IF( sec%ll_ice_section )THEN
1588
1589              IF ( transports_2d(1,jsec,jseg) .GE. 0.0 ) THEN
1590                 sec%transport(11,1) = sec%transport(11,1)+transports_2d(1,jsec,jseg)
1591              ELSE
1592                 sec%transport(12,1) = sec%transport(12,1)+transports_2d(1,jsec,jseg)
1593              ENDIF
1594
1595              IF ( transports_2d(3,jsec,jseg) .GE. 0.0 ) THEN
1596                 sec%transport(13,1) = sec%transport(13,1)+transports_2d(2,jsec,jseg)
1597              ELSE
1598                 sec%transport(14,1) = sec%transport(14,1)+transports_2d(2,jsec,jseg)
1599              ENDIF
1600
1601           ENDIF !end of ice case
1602#endif
1603 
1604        ENDDO !end of loop on the segment
1605
1606     ELSE  !if sec%nb_point =0
1607        sec%transport(1:2,:)=0.
1608        IF (sec%llstrpond) sec%transport(3:10,:)=0.
1609        IF (sec%ll_ice_section) sec%transport( 11:14,:)=0.
1610     ENDIF !end of sec%nb_point =0 case
1611
1612  END SUBROUTINE dia_dct_sum
1613 
1614  SUBROUTINE dia_dct_sum_h(sec,jsec)
1615     !!-------------------------------------------------------------
1616     !! Exactly as dia_dct_sum but for hourly files containing data summed at each time step
1617     !!
1618     !! Purpose: Average the transport over nn_dctwri time steps
1619     !! and sum over the density/salinity/temperature/depth classes
1620     !!
1621     !! Method:
1622     !!           Sum over relevant grid cells to obtain values
1623     !!           for each
1624     !!              There are several loops:                 
1625     !!              loop on the segment between 2 nodes
1626     !!              loop on the level jk
1627     !!              loop on the density/temperature/salinity/level classes
1628     !!              test on the density/temperature/salinity/level
1629     !!
1630     !!  ** Method  :Transport through a given section is equal to the sum of transports
1631     !!              computed on each proc.
1632     !!              On each proc,transport is equal to the sum of transport computed through
1633     !!              segments linking each point of sec%listPoint  with the next one.   
1634     !!
1635     !!-------------------------------------------------------------
1636     !! * arguments
1637     TYPE(SECTION),INTENT(INOUT) :: sec
1638     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section
1639
1640     TYPE(POINT_SECTION) :: k
1641     INTEGER  :: jk,jseg,jclass                        !loop on level/segment/classes
1642     REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point
1643     !!-------------------------------------------------------------
1644
1645     !! Sum the relevant segments to obtain values for each class
1646     IF(sec%nb_point .NE. 0)THEN   
1647
1648        !--------------------------------------!
1649        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  !
1650        !--------------------------------------!
1651        DO jseg=1,MAX(sec%nb_point-1,0)
1652           
1653           !-------------------------------------------------------------------------------------------
1654           ! Select the appropriate coordinate for computing the velocity of the segment
1655           !
1656           !                      CASE(0)                                    Case (2)
1657           !                      -------                                    --------
1658           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)     
1659           !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               |
1660           !                                                                            |
1661           !                                                                            |
1662           !                                                                            |
1663           !                      Case (3)                                            U(i,j)
1664           !                      --------                                              |
1665           !                                                                            |
1666           !  listPoint(jseg+1) F(i,j+1)                                                |
1667           !                        |                                                   |
1668           !                        |                                                   |
1669           !                        |                                 listPoint(jseg+1) F(i,j-1)
1670           !                        |                                           
1671           !                        |                                           
1672           !                     U(i,j+1)                                           
1673           !                        |                                       Case(1)     
1674           !                        |                                       ------     
1675           !                        |                                           
1676           !                        |                 listPoint(jseg+1)             listPoint(jseg)                           
1677           !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                           
1678           ! listPoint(jseg)     F(i,j)
1679           !
1680           !-------------------------------------------------------------------------------------------
1681
1682           SELECT CASE( sec%direction(jseg) )
1683           CASE(0)  ;   k = sec%listPoint(jseg)
1684           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
1685           CASE(2)  ;   k = sec%listPoint(jseg)
1686           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
1687           END SELECT
1688
1689           !---------------------------|
1690           !     LOOP ON THE LEVEL     |
1691           !---------------------------|
1692           !Sum of the transport on the vertical
1693           DO jk=1,mbathy(k%I,k%J)
1694
1695              ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
1696              SELECT CASE( sec%direction(jseg) )
1697              CASE(0,1)
1698                 ztn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) )
1699                 zsn   = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) )
1700                 zrhop = interp(k%I,k%J,jk,'V',rhop)
1701                 zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0)
1702                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1)
1703              CASE(2,3)
1704                 ztn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) )
1705                 zsn   = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) )
1706                 zrhop = interp(k%I,k%J,jk,'U',rhop)
1707                 zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0)
1708                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1) 
1709              END SELECT
1710
1711              zfsdep= fsdept(k%I,k%J,jk)
1712 
1713              !-------------------------------
1714              !  LOOP ON THE DENSITY CLASSES |
1715              !-------------------------------
1716              !The computation is made for each density/heat/salt/... class
1717              DO jclass=1,MAX(1,sec%nb_class-1)
1718
1719                 !----------------------------------------------!
1720                 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL!
1721                 !----------------------------------------------!
1722 
1723                 IF ( (                                                    &
1724                    ((( zrhop .GE. (sec%zsigp(jclass)+1000.  )) .AND.      &
1725                    (   zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR.     &
1726                    ( sec%zsigp(jclass) .EQ. 99.)) .AND.                   &
1727
1728                    ((( zrhoi .GE. (sec%zsigi(jclass) + 1000.  )) .AND.    &
1729                    (   zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR.     &
1730                    ( sec%zsigi(jclass) .EQ. 99.)) .AND.                   &
1731
1732                    ((( zsn .GT. sec%zsal(jclass)) .AND.                   &
1733                    (   zsn .LE. sec%zsal(jclass+1))) .OR.                 &
1734                    ( sec%zsal(jclass) .EQ. 99.)) .AND.                    &
1735
1736                    ((( ztn .GE. sec%ztem(jclass)) .AND.                   &
1737                    (   ztn .LE. sec%ztem(jclass+1))) .OR.                 &
1738                    ( sec%ztem(jclass) .EQ.99.)) .AND.                     &
1739
1740                    ((( zfsdep .GE. sec%zlay(jclass)) .AND.                &
1741                    (   zfsdep .LE. sec%zlay(jclass+1))) .OR.              &
1742                    ( sec%zlay(jclass) .EQ. 99. ))                         &
1743                                                                   ))   THEN
1744
1745                    !SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS
1746                    !----------------------------------------------------------------------------
1747                    IF (transports_3d_h(1,jsec,jseg,jk) .GE. 0.0) THEN
1748                       sec%transport_h(1,jclass) = sec%transport_h(1,jclass)+transports_3d_h(1,jsec,jseg,jk)
1749                    ELSE
1750                       sec%transport_h(2,jclass) = sec%transport_h(2,jclass)+transports_3d_h(1,jsec,jseg,jk)
1751                    ENDIF
1752                    IF( sec%llstrpond )THEN
1753
1754                       IF( transports_3d_h(1,jsec,jseg,jk) .NE. 0._wp ) THEN
1755
1756                          IF (transports_3d_h(2,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) .GE. 0.0 ) THEN
1757                             sec%transport_h(3,jclass) = sec%transport_h(3,jclass)+transports_3d_h(2,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk)
1758                          ELSE
1759                             sec%transport_h(4,jclass) = sec%transport_h(4,jclass)+transports_3d_h(2,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk)
1760                          ENDIF
1761
1762                          IF ( transports_3d_h(3,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) .GE. 0.0 ) THEN
1763                             sec%transport_h(5,jclass) = sec%transport_h(5,jclass)+transports_3d_h(3,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk)
1764                          ELSE
1765                             sec%transport_h(6,jclass) = sec%transport_h(6,jclass)+transports_3d_h(3,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk)
1766                          ENDIF
1767
1768                       ENDIF
1769
1770                       IF ( transports_3d_h(4,jsec,jseg,jk) .GE. 0.0 ) THEN
1771                          sec%transport_h(7,jclass) = sec%transport_h(7,jclass)+transports_3d_h(4,jsec,jseg,jk)
1772                       ELSE
1773                          sec%transport_h(8,jclass) = sec%transport_h(8,jclass)+transports_3d_h(4,jsec,jseg,jk)
1774                       ENDIF
1775
1776                       IF ( transports_3d_h(5,jsec,jseg,jk) .GE. 0.0 ) THEN
1777                          sec%transport_h( 9,jclass) = sec%transport_h( 9,jclass)+transports_3d_h(5,jsec,jseg,jk)
1778                       ELSE
1779                          sec%transport_h(10,jclass) = sec%transport_h(10,jclass)+transports_3d_h(5,jsec,jseg,jk)
1780                       ENDIF
1781
1782                    ELSE
1783                       sec%transport_h( 3,jclass) = 0._wp
1784                       sec%transport_h( 4,jclass) = 0._wp
1785                       sec%transport_h( 5,jclass) = 0._wp
1786                       sec%transport_h( 6,jclass) = 0._wp
1787                       sec%transport_h( 7,jclass) = 0._wp
1788                       sec%transport_h( 8,jclass) = 0._wp
1789                       sec%transport_h( 9,jclass) = 0._wp
1790                       sec%transport_h(10,jclass) = 0._wp
1791                    ENDIF
1792
1793                 ENDIF ! end of test if point is in class
1794   
1795              ENDDO ! end of loop on the classes
1796
1797           ENDDO ! loop over jk
1798
1799#if defined key_lim2 || defined key_lim3
1800
1801           !ICE CASE   
1802           IF( sec%ll_ice_section )THEN
1803
1804              IF ( transports_2d_h(1,jsec,jseg) .GE. 0.0 ) THEN
1805                 sec%transport_h(11,1) = sec%transport_h(11,1)+transports_2d_h(1,jsec,jseg)
1806              ELSE
1807                 sec%transport_h(12,1) = sec%transport_h(12,1)+transports_2d_h(1,jsec,jseg)
1808              ENDIF
1809
1810              IF ( transports_2d_h(3,jsec,jseg) .GE. 0.0 ) THEN
1811                 sec%transport_h(13,1) = sec%transport_h(13,1)+transports_2d_h(2,jsec,jseg)
1812              ELSE
1813                 sec%transport_h(14,1) = sec%transport_h(14,1)+transports_2d_h(2,jsec,jseg)
1814              ENDIF
1815
1816           ENDIF !end of ice case
1817#endif
1818 
1819        ENDDO !end of loop on the segment
1820
1821     ELSE  !if sec%nb_point =0
1822        sec%transport_h(1:2,:)=0.
1823        IF (sec%llstrpond) sec%transport_h(3:10,:)=0.
1824        IF (sec%ll_ice_section) sec%transport_h( 11:14,:)=0.
1825     ENDIF !end of sec%nb_point =0 case
1826
1827  END SUBROUTINE dia_dct_sum_h
1828 
1829  SUBROUTINE dia_dct_wri_NOOS(kt,ksec,sec)
1830     !!-------------------------------------------------------------
1831     !! Write transport output in numdct using NOOS formatting
1832     !!
1833     !! Purpose: Write  transports in ascii files
1834     !!
1835     !! Method:
1836     !!        1. Write volume transports in "volume_transport"
1837     !!           Unit: Sv : area * Velocity / 1.e6
1838     !!
1839     !!        2. Write heat transports in "heat_transport"
1840     !!           Unit: Peta W : area * Velocity * T * rhau * Cp / 1.e15
1841     !!
1842     !!        3. Write salt transports in "salt_transport"
1843     !!           Unit: 10^9 g m^3 / s : area * Velocity * S / 1.e6
1844     !!
1845     !!-------------------------------------------------------------
1846     !!arguments
1847     INTEGER, INTENT(IN)          :: kt          ! time-step
1848     TYPE(SECTION), INTENT(INOUT) :: sec         ! section to write   
1849     INTEGER ,INTENT(IN)          :: ksec        ! section number
1850
1851     !!local declarations
1852     INTEGER               :: jclass,ji             ! Dummy loop
1853     CHARACTER(len=2)      :: classe             ! Classname
1854     
1855     REAL(wp)              :: zbnd1,zbnd2        ! Class bounds
1856     REAL(wp)              :: zslope             ! section's slope coeff
1857     !
1858     REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace
1859     CHARACTER(len=3)      :: noos_sect_name             ! Classname
1860     CHARACTER(len=25)      :: noos_var_sect_name
1861     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   noos_iom_dummy
1862     INTEGER               :: IERR
1863     
1864     REAL(wp), DIMENSION(3) :: tmp_iom_output
1865     REAL(wp)               :: max_iom_val
1866     
1867     !!-------------------------------------------------------------
1868     
1869     
1870     
1871     IF( lwp ) THEN
1872        WRITE(numout,*) " "
1873        WRITE(numout,*) "dia_dct_wri_NOOS: write transports through sections at timestep: ", kt
1874        WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
1875     ENDIF   
1876       
1877     CALL wrk_alloc(nb_type , zsumclasses ) 
1878
1879     zsumclasses(:)=0._wp
1880     zslope = sec%slopeSection       
1881     
1882     IF( lwp ) THEN
1883         IF  ( ln_dct_ascii ) THEN
1884             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
1885         ELSE
1886             WRITE(numdct_NOOS) nyear,nmonth,nday,ksec-1,sec%nb_class-1,sec%name
1887         ENDIF 
1888     ENDIF
1889   
1890     ! Sum all classes together, to give one values per type (pos tran, neg vol trans etc...).
1891     DO jclass=1,MAX(1,sec%nb_class-1)
1892        zsumclasses(1:nb_type)=zsumclasses(1:nb_type)+sec%transport(1:nb_type,jclass)
1893     ENDDO
1894 
1895     classe   = 'total   '
1896     zbnd1   = 0._wp
1897     zbnd2   = 0._wp
1898     
1899     
1900     
1901     write (noos_sect_name, "(I0.3)")  ksec
1902     
1903     IF ( ln_dct_iom_cont ) THEN
1904         max_iom_val = 1.e10
1905         ALLOCATE( noos_iom_dummy(jpi,jpj,3),  STAT= ierr )
1906            IF( ierr /= 0 )   CALL ctl_stop( 'dia_dct_wri_NOOS: failed to allocate noos_iom_dummy array' )
1907     ENDIF
1908     
1909!     JT
1910!     JT
1911!     JT
1912!     JT   I think changing the sign on the output based on the zslope value is redunant.
1913!     JT
1914!     JT
1915!     JT
1916!     JT
1917!     JT
1918!
1919!
1920!
1921!     IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN
1922!         
1923!         IF( lwp ) THEN
1924!             WRITE(numdct_NOOS,'(9e12.4E2)') -(zsumclasses( 1)+zsumclasses( 2)), -zsumclasses( 2),-zsumclasses( 1),   &
1925!                                            -(zsumclasses( 7)+zsumclasses( 8)), -zsumclasses( 8),-zsumclasses( 7),   &
1926!                                            -(zsumclasses( 9)+zsumclasses(10)), -zsumclasses(10),-zsumclasses( 9)
1927!             CALL FLUSH(numdct_NOOS)
1928!         endif
1929
1930!         
1931!         IF ( ln_dct_iom_cont ) THEN
1932!         
1933!             noos_iom_dummy(:,:,:) = 0.
1934!             
1935!             tmp_iom_output(:) = 0.
1936!             tmp_iom_output(1) = -(zsumclasses( 1)+zsumclasses( 2))
1937!             tmp_iom_output(2) = -zsumclasses( 2)
1938!             tmp_iom_output(3) = -zsumclasses( 1)
1939!             
1940!             ! Convert to Sv
1941!             tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
1942!             tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
1943!             tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
1944!             
1945!             ! limit maximum and minimum values in iom_put
1946!             if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
1947!             if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
1948!             if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
1949!             if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
1950!             if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
1951!             if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
1952!             
1953!             ! Set NaN's to Zero         
1954!             if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
1955!             if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
1956!             if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
1957!             
1958!             noos_iom_dummy(:,:,1) = tmp_iom_output(1)
1959!             noos_iom_dummy(:,:,2) = tmp_iom_output(2)
1960!             noos_iom_dummy(:,:,3) = tmp_iom_output(3)
1961!             
1962!             noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_trans'
1963!             if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
1964!             CALL iom_put( noos_var_sect_name,  noos_iom_dummy )
1965!             noos_iom_dummy(:,:,:) = 0.         
1966!             tmp_iom_output(:) = 0.
1967!             tmp_iom_output(1) = -(zsumclasses( 7)+zsumclasses( 8))
1968!             tmp_iom_output(2) = -zsumclasses( 8)
1969!             tmp_iom_output(3) = -zsumclasses( 7)
1970!             
1971!             ! Convert to TJ/s
1972!             tmp_iom_output(1) = tmp_iom_output(1)*1.E-12
1973!             tmp_iom_output(2) = tmp_iom_output(2)*1.E-12
1974!             tmp_iom_output(3) = tmp_iom_output(3)*1.E-12
1975!             
1976!             ! limit maximum and minimum values in iom_put
1977!             if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
1978!             if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
1979!             if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
1980!             if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
1981!             if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
1982!             if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
1983!             
1984!             ! Set NaN's to Zero         
1985!             if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
1986!             if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
1987!             if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
1988!             
1989!             noos_iom_dummy(:,:,1) = tmp_iom_output(1)
1990!             noos_iom_dummy(:,:,2) = tmp_iom_output(2)
1991!             noos_iom_dummy(:,:,3) = tmp_iom_output(3)
1992!             
1993!             !noos_iom_dummy(:,:,1) = -(zsumclasses( 7)+zsumclasses( 8))
1994!             !noos_iom_dummy(:,:,2) = -zsumclasses( 8)
1995!             !noos_iom_dummy(:,:,3) = -zsumclasses( 7)
1996!             
1997!             noos_var_sect_name =  "noos_" // trim(noos_sect_name) // '_heat'
1998!             if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
1999!             CALL iom_put( noos_var_sect_name,  noos_iom_dummy )
2000!             
2001!             noos_iom_dummy(:,:,:) = 0.         
2002!             tmp_iom_output(:) = 0.         
2003!             tmp_iom_output(1) = -(zsumclasses( 9)+zsumclasses( 10))
2004!             tmp_iom_output(2) = -zsumclasses( 10)
2005!             tmp_iom_output(3) = -zsumclasses( 9)
2006!             
2007!             ! Convert to  MT/s
2008!             tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
2009!             tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
2010!             tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
2011!             
2012!             ! limit maximum and minimum values in iom_put
2013!             if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
2014!             if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
2015!             if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
2016!             if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
2017!             if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
2018!             if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
2019!             
2020!             ! Set NaN's to Zero         
2021!             if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
2022!             if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
2023!             if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
2024!             
2025!             noos_iom_dummy(:,:,1) = tmp_iom_output(1)
2026!             noos_iom_dummy(:,:,2) = tmp_iom_output(2)
2027!             noos_iom_dummy(:,:,3) = tmp_iom_output(3)
2028!             
2029!             !noos_iom_dummy(:,:,1) = -(zsumclasses( 9)+zsumclasses( 10))
2030!             !noos_iom_dummy(:,:,2) = -zsumclasses( 10)
2031!             !noos_iom_dummy(:,:,3) = -zsumclasses( 9)
2032!             
2033!             noos_var_sect_name =  "noos_" // trim(noos_sect_name) // '_salt'
2034!             if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
2035!             CALL iom_put( noos_var_sect_name,  noos_iom_dummy )
2036!             noos_iom_dummy(:,:,:) = 0.     
2037!             tmp_iom_output(:) = 0.   
2038!         ENDIF
2039!     ELSE
2040!        IF( lwp ) THEN
2041!            WRITE(numdct_NOOS,'(9e12.4E2)')   zsumclasses( 1)+zsumclasses( 2) ,  zsumclasses( 1), zsumclasses( 2),   &
2042!                                            zsumclasses( 7)+zsumclasses( 8) ,  zsumclasses( 7), zsumclasses( 8),   &
2043!                                            zsumclasses( 9)+zsumclasses(10) ,  zsumclasses( 9), zsumclasses(10)
2044!            CALL FLUSH(numdct_NOOS)
2045!        endif
2046!         
2047!         
2048!        IF ( ln_dct_iom_cont ) THEN
2049!           
2050!             noos_iom_dummy(:,:,:) = 0.
2051!             tmp_iom_output(:) = 0.
2052!             
2053!             tmp_iom_output(1) = (zsumclasses( 1)+zsumclasses( 2))
2054!             tmp_iom_output(2) = zsumclasses( 1)
2055!             tmp_iom_output(3) = zsumclasses( 2)
2056!             
2057!             ! Convert to Sv
2058!             tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
2059!             tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
2060!             tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
2061!             
2062!             ! limit maximum and minimum values in iom_put
2063!             if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
2064!             if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
2065!             if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
2066!             if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
2067!             if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
2068!             if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
2069!             
2070!             ! Set NaN's to Zero         
2071!             if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
2072!             if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
2073!             if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
2074!             
2075!             noos_iom_dummy(:,:,1) = tmp_iom_output(1)
2076!             noos_iom_dummy(:,:,2) = tmp_iom_output(2)
2077!             noos_iom_dummy(:,:,3) = tmp_iom_output(3)
2078!             
2079!             !noos_iom_dummy(:,:,1) = (zsumclasses( 1)+zsumclasses( 2))
2080!             !noos_iom_dummy(:,:,2) = zsumclasses( 1)
2081!             !noos_iom_dummy(:,:,3) = zsumclasses( 2)
2082!             
2083!             
2084!             
2085!             noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_trans'
2086!             if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
2087!             CALL iom_put( noos_var_sect_name,  noos_iom_dummy )
2088!             noos_iom_dummy(:,:,:) = 0.
2089!             tmp_iom_output(:) = 0.
2090!             
2091!             tmp_iom_output(1) = (zsumclasses( 7)+zsumclasses( 8))
2092!             tmp_iom_output(2) = zsumclasses( 7)
2093!             tmp_iom_output(3) = zsumclasses( 8)
2094!             
2095!             ! Convert to TJ/s
2096!             tmp_iom_output(1) = tmp_iom_output(1)*1.E-12
2097!             tmp_iom_output(2) = tmp_iom_output(2)*1.E-12
2098!             tmp_iom_output(3) = tmp_iom_output(3)*1.E-12
2099!             
2100!             ! limit maximum and minimum values in iom_put
2101!             if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
2102!             if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
2103!             if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
2104!             if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
2105!             if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
2106!             if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
2107!             
2108!             ! Set NaN's to Zero         
2109!             if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
2110!             if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
2111!             if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
2112!             
2113!             noos_iom_dummy(:,:,1) = tmp_iom_output(1)
2114!             noos_iom_dummy(:,:,2) = tmp_iom_output(2)
2115!             noos_iom_dummy(:,:,3) = tmp_iom_output(3)
2116!             
2117!             !noos_iom_dummy(:,:,1) = (zsumclasses( 7)+zsumclasses( 8))
2118!             !noos_iom_dummy(:,:,2) = zsumclasses( 7)
2119!             !noos_iom_dummy(:,:,3) = zsumclasses( 8)
2120!             
2121!             noos_var_sect_name =  "noos_" // trim(noos_sect_name) // '_heat'
2122!             if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
2123!             CALL iom_put(noos_var_sect_name,  noos_iom_dummy ) 
2124!             noos_iom_dummy(:,:,:) = 0.
2125!             tmp_iom_output(:) = 0.
2126!             
2127!             tmp_iom_output(1) = (zsumclasses( 9)+zsumclasses( 10))
2128!             tmp_iom_output(2) = zsumclasses( 9)
2129!             tmp_iom_output(3) = zsumclasses( 10)
2130!             
2131!             ! Convert to  MT/s
2132!             tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
2133!             tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
2134!             tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
2135!             
2136!             
2137!             ! limit maximum and minimum values in iom_put
2138!             if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
2139!             if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
2140!             if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
2141!             if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
2142!             if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
2143!             if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
2144!             
2145!             ! Set NaN's to Zero         
2146!             if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
2147!             if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
2148!             if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
2149!             
2150!             noos_iom_dummy(:,:,1) = tmp_iom_output(1)
2151!             noos_iom_dummy(:,:,2) = tmp_iom_output(2)
2152!             noos_iom_dummy(:,:,3) = tmp_iom_output(3)
2153!             
2154!             !noos_iom_dummy(:,:,1) = (zsumclasses( 9)+zsumclasses( 10))
2155!             !noos_iom_dummy(:,:,2) = zsumclasses( 9)
2156!             !noos_iom_dummy(:,:,3) = zsumclasses( 10)
2157!             
2158!             noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_salt'
2159!             if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
2160!             CALL iom_put(noos_var_sect_name,  noos_iom_dummy )
2161!             noos_iom_dummy(:,:,:) = 0.         
2162!             tmp_iom_output(:) = 0.
2163!         ENDIF
2164!         
2165!     ENDIF
2166     
2167
2168
2169
2170
2171
2172
2173
2174     
2175     
2176     IF( lwp ) THEN
2177        IF  ( ln_dct_ascii ) THEN
2178             !WRITE(numdct_NOOS,'(9e12.4E2)')   zsumclasses( 1)+zsumclasses( 2) ,  zsumclasses( 1), zsumclasses( 2),   &
2179             WRITE(numdct_NOOS,'(3F18.3,6e16.8E2)')   zsumclasses( 1)+zsumclasses( 2) ,  zsumclasses( 1), zsumclasses( 2),   &
2180                                             zsumclasses( 7)+zsumclasses( 8) ,  zsumclasses( 7), zsumclasses( 8),   &
2181                                             zsumclasses( 9)+zsumclasses(10) ,  zsumclasses( 9), zsumclasses(10)
2182             CALL FLUSH(numdct_NOOS)
2183        ELSE
2184             WRITE(numdct_NOOS)   zsumclasses( 1)+zsumclasses( 2) ,  zsumclasses( 1), zsumclasses( 2),   &
2185                                  zsumclasses( 7)+zsumclasses( 8) ,  zsumclasses( 7), zsumclasses( 8),   &
2186                                  zsumclasses( 9)+zsumclasses(10) ,  zsumclasses( 9), zsumclasses(10)
2187             CALL FLUSH(numdct_NOOS) 
2188         ENDIF
2189     ENDIF
2190         
2191    IF ( ln_dct_iom_cont ) THEN
2192       
2193         noos_iom_dummy(:,:,:) = 0.
2194         tmp_iom_output(:) = 0.
2195         
2196         tmp_iom_output(1) = (zsumclasses( 1)+zsumclasses( 2))
2197         tmp_iom_output(2) = zsumclasses( 1)
2198         tmp_iom_output(3) = zsumclasses( 2)
2199         
2200         ! Convert to Sv
2201         tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
2202         tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
2203         tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
2204         
2205         ! limit maximum and minimum values in iom_put
2206         if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
2207         if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
2208         if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
2209         if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
2210         if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
2211         if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
2212         
2213         ! Set NaN's to Zero         
2214         if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
2215         if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
2216         if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
2217         
2218         noos_iom_dummy(:,:,1) = tmp_iom_output(1)
2219         noos_iom_dummy(:,:,2) = tmp_iom_output(2)
2220         noos_iom_dummy(:,:,3) = tmp_iom_output(3)
2221         
2222         !noos_iom_dummy(:,:,1) = (zsumclasses( 1)+zsumclasses( 2))
2223         !noos_iom_dummy(:,:,2) = zsumclasses( 1)
2224         !noos_iom_dummy(:,:,3) = zsumclasses( 2)
2225         
2226         
2227         
2228         noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_trans'
2229         if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
2230         CALL iom_put( noos_var_sect_name,  noos_iom_dummy )
2231         noos_iom_dummy(:,:,:) = 0.
2232         tmp_iom_output(:) = 0.
2233         
2234         tmp_iom_output(1) = (zsumclasses( 7)+zsumclasses( 8))
2235         tmp_iom_output(2) = zsumclasses( 7)
2236         tmp_iom_output(3) = zsumclasses( 8)
2237         
2238         ! Convert to TJ/s
2239         tmp_iom_output(1) = tmp_iom_output(1)*1.E-12
2240         tmp_iom_output(2) = tmp_iom_output(2)*1.E-12
2241         tmp_iom_output(3) = tmp_iom_output(3)*1.E-12
2242         
2243         ! limit maximum and minimum values in iom_put
2244         if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
2245         if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
2246         if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
2247         if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
2248         if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
2249         if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
2250         
2251         ! Set NaN's to Zero         
2252         if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
2253         if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
2254         if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
2255         
2256         noos_iom_dummy(:,:,1) = tmp_iom_output(1)
2257         noos_iom_dummy(:,:,2) = tmp_iom_output(2)
2258         noos_iom_dummy(:,:,3) = tmp_iom_output(3)
2259         
2260         !noos_iom_dummy(:,:,1) = (zsumclasses( 7)+zsumclasses( 8))
2261         !noos_iom_dummy(:,:,2) = zsumclasses( 7)
2262         !noos_iom_dummy(:,:,3) = zsumclasses( 8)
2263         
2264         noos_var_sect_name =  "noos_" // trim(noos_sect_name) // '_heat'
2265         if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
2266         CALL iom_put(noos_var_sect_name,  noos_iom_dummy ) 
2267         noos_iom_dummy(:,:,:) = 0.
2268         tmp_iom_output(:) = 0.
2269         
2270         tmp_iom_output(1) = (zsumclasses( 9)+zsumclasses( 10))
2271         tmp_iom_output(2) = zsumclasses( 9)
2272         tmp_iom_output(3) = zsumclasses( 10)
2273         
2274         ! Convert to  MT/s
2275         tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
2276         tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
2277         tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
2278         
2279         
2280         ! limit maximum and minimum values in iom_put
2281         if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
2282         if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
2283         if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
2284         if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
2285         if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
2286         if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
2287         
2288         ! Set NaN's to Zero         
2289         if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
2290         if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
2291         if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
2292         
2293         noos_iom_dummy(:,:,1) = tmp_iom_output(1)
2294         noos_iom_dummy(:,:,2) = tmp_iom_output(2)
2295         noos_iom_dummy(:,:,3) = tmp_iom_output(3)
2296         
2297         !noos_iom_dummy(:,:,1) = (zsumclasses( 9)+zsumclasses( 10))
2298         !noos_iom_dummy(:,:,2) = zsumclasses( 9)
2299         !noos_iom_dummy(:,:,3) = zsumclasses( 10)
2300         
2301         noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_salt'
2302         if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
2303         CALL iom_put(noos_var_sect_name,  noos_iom_dummy )
2304         noos_iom_dummy(:,:,:) = 0.         
2305         tmp_iom_output(:) = 0.
2306
2307
2308        DEALLOCATE(noos_iom_dummy)
2309     ENDIF
2310     
2311
2312     DO jclass=1,MAX(1,sec%nb_class-1)
2313   
2314        classe   = 'N       '
2315        zbnd1   = 0._wp
2316        zbnd2   = 0._wp
2317
2318        !insitu density classes transports
2319        IF( ( sec%zsigi(jclass)   .NE. 99._wp ) .AND. &
2320            ( sec%zsigi(jclass+1) .NE. 99._wp )       )THEN
2321           classe = 'DI       '
2322           zbnd1 = sec%zsigi(jclass)
2323           zbnd2 = sec%zsigi(jclass+1)
2324        ENDIF
2325        !potential density classes transports
2326        IF( ( sec%zsigp(jclass)   .NE. 99._wp ) .AND. &
2327            ( sec%zsigp(jclass+1) .NE. 99._wp )       )THEN
2328           classe = 'DP      '
2329           zbnd1 = sec%zsigp(jclass)
2330           zbnd2 = sec%zsigp(jclass+1)
2331        ENDIF
2332        !depth classes transports
2333        IF( ( sec%zlay(jclass)    .NE. 99._wp ) .AND. &
2334            ( sec%zlay(jclass+1)  .NE. 99._wp )       )THEN
2335           classe = 'Z       '
2336           zbnd1 = sec%zlay(jclass)
2337           zbnd2 = sec%zlay(jclass+1)
2338        ENDIF
2339        !salinity classes transports
2340        IF( ( sec%zsal(jclass) .NE. 99._wp    ) .AND. &
2341            ( sec%zsal(jclass+1) .NE. 99._wp  )       )THEN
2342           classe = 'S       '
2343           zbnd1 = sec%zsal(jclass)
2344           zbnd2 = sec%zsal(jclass+1)   
2345        ENDIF
2346        !temperature classes transports
2347        IF( ( sec%ztem(jclass) .NE. 99._wp     ) .AND. &
2348            ( sec%ztem(jclass+1) .NE. 99._wp     )       ) THEN
2349           classe = 'T       '
2350           zbnd1 = sec%ztem(jclass)
2351           zbnd2 = sec%ztem(jclass+1)
2352        ENDIF
2353                 
2354        !write volume transport per class
2355        IF( lwp ) THEN
2356           
2357            IF  ( ln_dct_ascii ) THEN
2358                CALL FLUSH(numdct_NOOS) ! JT crash
2359
2360
2361
2362                !JT    IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN
2363                !JT      WRITE(numdct_NOOS,'(9e12.4E2)') -(sec%transport( 1,jclass)+sec%transport( 2,jclass)),-sec%transport( 2,jclass),-sec%transport( 1,jclass), &
2364                !JT                                      -(sec%transport( 7,jclass)+sec%transport( 8,jclass)),-sec%transport( 8,jclass),-sec%transport( 7,jclass), &
2365                !JT                                      -(sec%transport( 9,jclass)+sec%transport(10,jclass)),-sec%transport(10,jclass),-sec%transport( 9,jclass)
2366                !JT    ELSE
2367                !JT      WRITE(numdct_NOOS,'(9e12.4E2)')   sec%transport( 1,jclass)+sec%transport( 2,jclass) , sec%transport( 1,jclass), sec%transport( 2,jclass), &
2368                !JT                                        sec%transport( 7,jclass)+sec%transport( 8,jclass) , sec%transport( 7,jclass), sec%transport( 8,jclass), &
2369                !JT                                        sec%transport( 9,jclass)+sec%transport(10,jclass) , sec%transport( 9,jclass), sec%transport(10,jclass)
2370                !JT    ENDIF
2371
2372
2373                !WRITE(numdct_NOOS,'(9e12.4E2)')   sec%transport( 1,jclass)+sec%transport( 2,jclass) , sec%transport( 1,jclass), sec%transport( 2,jclass), &
2374                !                                  sec%transport( 7,jclass)+sec%transport( 8,jclass) , sec%transport( 7,jclass), sec%transport( 8,jclass), &
2375                !                                  sec%transport( 9,jclass)+sec%transport(10,jclass) , sec%transport( 9,jclass), sec%transport(10,jclass)
2376                WRITE(numdct_NOOS,'(3F18.3,6e16.8E2)')   sec%transport( 1,jclass)+sec%transport( 2,jclass) , sec%transport( 1,jclass), sec%transport( 2,jclass), &
2377                                                         sec%transport( 7,jclass)+sec%transport( 8,jclass) , sec%transport( 7,jclass), sec%transport( 8,jclass), &
2378                                                         sec%transport( 9,jclass)+sec%transport(10,jclass) , sec%transport( 9,jclass), sec%transport(10,jclass)
2379            ELSE
2380
2381                CALL FLUSH(numdct_NOOS) ! JT crash
2382                WRITE(numdct_NOOS)   sec%transport( 1,jclass)+sec%transport( 2,jclass) , sec%transport( 1,jclass), sec%transport( 2,jclass), &
2383                                     sec%transport( 7,jclass)+sec%transport( 8,jclass) , sec%transport( 7,jclass), sec%transport( 8,jclass), &
2384                                     sec%transport( 9,jclass)+sec%transport(10,jclass) , sec%transport( 9,jclass), sec%transport(10,jclass)
2385            ENDIF
2386        ENDIF
2387
2388     ENDDO
2389     
2390     !IF  ( ln_dct_ascii ) THEN
2391        if ( lwp ) CALL FLUSH(numdct_NOOS)
2392     !ENDIF
2393
2394     CALL wrk_dealloc(nb_type , zsumclasses ) 
2395
2396  END SUBROUTINE dia_dct_wri_NOOS
2397
2398  SUBROUTINE dia_dct_wri_NOOS_h(hr,ksec,sec)
2399     !!-------------------------------------------------------------
2400     !! As routine dia_dct_wri_NOOS but for hourly output files
2401     !!
2402     !! Write transport output in numdct using NOOS formatting
2403     !!
2404     !! Purpose: Write  transports in ascii files
2405     !!
2406     !! Method:
2407     !!        1. Write volume transports in "volume_transport"
2408     !!           Unit: Sv : area * Velocity / 1.e6
2409     !!
2410     !!-------------------------------------------------------------
2411     !!arguments
2412     INTEGER, INTENT(IN)          :: hr          ! hour => effectively kt/12
2413     TYPE(SECTION), INTENT(INOUT) :: sec         ! section to write   
2414     INTEGER ,INTENT(IN)          :: ksec        ! section number
2415
2416     !!local declarations
2417     INTEGER               :: jclass,jhr            ! Dummy loop
2418     CHARACTER(len=2)      :: classe             ! Classname
2419     REAL(wp)              :: zbnd1,zbnd2        ! Class bounds
2420     REAL(wp)              :: zslope             ! section's slope coeff
2421     !
2422     REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace
2423     CHARACTER(len=3)      :: noos_sect_name             ! Classname
2424     CHARACTER(len=25)      :: noos_var_sect_name
2425     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   noos_iom_dummy
2426     INTEGER               :: IERR
2427     
2428     !!-------------------------------------------------------------
2429
2430     IF( lwp ) THEN
2431        WRITE(numout,*) " "
2432        WRITE(numout,*) "dia_dct_wri_NOOS_h: write transports through section Transect:",ksec-1," at timestep: ", hr
2433        WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
2434     ENDIF   
2435     
2436     CALL wrk_alloc(nb_type , zsumclasses ) 
2437     
2438     
2439     write (noos_sect_name, "(I03)")  ksec
2440     
2441     ALLOCATE( noos_iom_dummy(jpi,jpj,3),  STAT= ierr )
2442        IF( ierr /= 0 )   CALL ctl_stop( 'dia_dct_wri_NOOS_h: failed to allocate noos_iom_dummy array' )
2443
2444
2445
2446     
2447
2448     zsumclasses(:)=0._wp
2449     zslope = sec%slopeSection       
2450
2451     ! Sum up all classes, to give the total per type (pos vol trans, neg vol trans etc...)
2452     DO jclass=1,MAX(1,sec%nb_class-1)
2453        zsumclasses(1:nb_type)=zsumclasses(1:nb_type)+sec%transport_h(1:nb_type,jclass)
2454     ENDDO
2455       
2456     
2457     ! JT
2458     ! JT
2459     ! JT
2460     ! JT
2461     ! JT I think changing the sign of output according to the zslope is redundant
2462     ! JT
2463     ! JT
2464     ! JT
2465     ! JT
2466
2467     ! JT     !write volume transport per class
2468     ! JT     ! Sum positive and vol trans for all classes in first cell of array
2469     ! JT     IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN
2470     ! JT        ! JT z_hr_output(ksec,hr,1)=-(zsumclasses(1)+zsumclasses(2))
2471     ! JT        z_hr_output(ksec,1,1)=-(zsumclasses(1)+zsumclasses(2))
2472     ! JT     ELSE
2473     ! JT        ! JT z_hr_output(ksec,hr,1)= (zsumclasses(1)+zsumclasses(2))
2474     ! JT        z_hr_output(ksec,1,1)= (zsumclasses(1)+zsumclasses(2))
2475     ! JT     ENDIF
2476     ! JT
2477     ! JT     ! Sum positive and vol trans for each classes in following cell of array
2478     ! JT     DO jclass=1,MAX(1,sec%nb_class-1)
2479     ! JT        IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN
2480     ! JT           ! JT z_hr_output(ksec,hr,jclass+1)=-(sec%transport_h(1,jclass)+sec%transport_h(2,jclass))
2481     ! JT           z_hr_output(ksec,1,jclass+1)=-(sec%transport_h(1,jclass)+sec%transport_h(2,jclass))
2482     ! JT        ELSE
2483     ! JT           ! JT z_hr_output(ksec,hr,jclass+1)= (sec%transport_h(1,jclass)+sec%transport_h(2,jclass))
2484     ! JT           z_hr_output(ksec,1,jclass+1)= (sec%transport_h(1,jclass)+sec%transport_h(2,jclass))
2485     ! JT        ENDIF
2486     ! JT     ENDDO
2487     
2488     !write volume transport per class
2489     ! Sum positive and vol trans for all classes in first cell of array
2490
2491     z_hr_output(ksec,1,1)= (zsumclasses(1)+zsumclasses(2))
2492     z_hr_output(ksec,2,1)= zsumclasses(1)
2493     z_hr_output(ksec,3,1)= zsumclasses(2)
2494
2495     ! Sum positive and vol trans for each classes in following cell of array
2496     DO jclass=1,MAX(1,sec%nb_class-1)
2497        z_hr_output(ksec,1,jclass+1)= (sec%transport_h(1,jclass)+sec%transport_h(2,jclass))
2498        z_hr_output(ksec,2,jclass+1)= sec%transport_h(1,jclass)
2499        z_hr_output(ksec,3,jclass+1)= sec%transport_h(2,jclass)
2500     ENDDO
2501
2502   
2503    IF( lwp )  THEN
2504        ! JT IF ( hr .eq. 48._wp ) THEN
2505        ! 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
2506        ! JT    DO jhr=25,48
2507        ! 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) )
2508        ! JT    ENDDO
2509        ! JT ENDIF
2510
2511
2512
2513        IF ( ln_dct_ascii ) THEN
2514            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
2515            WRITE(numdct_NOOS_h,'(11F18.3)')  z_hr_output(ksec,1,1), (z_hr_output(ksec,1,jclass+1), jclass=1,MAX(1,10) )
2516            WRITE(numdct_NOOS_h,'(11F18.3)')  z_hr_output(ksec,2,1), (z_hr_output(ksec,2,jclass+1), jclass=1,MAX(1,10) )
2517            WRITE(numdct_NOOS_h,'(11F18.3)')  z_hr_output(ksec,3,1), (z_hr_output(ksec,3,jclass+1), jclass=1,MAX(1,10) )
2518            CALL FLUSH(numdct_NOOS_h)
2519        ELSE
2520            WRITE(numdct_NOOS_h) nyear,nmonth,nday,MOD(hr,24),ksec-1,sec%nb_class-1
2521            WRITE(numdct_NOOS_h)  z_hr_output(ksec,1,1), (z_hr_output(ksec,1,jclass+1), jclass=1,MAX(1,10) )
2522            WRITE(numdct_NOOS_h)  z_hr_output(ksec,2,1), (z_hr_output(ksec,2,jclass+1), jclass=1,MAX(1,10) )
2523            WRITE(numdct_NOOS_h)  z_hr_output(ksec,3,1), (z_hr_output(ksec,3,jclass+1), jclass=1,MAX(1,10) )
2524            CALL FLUSH(numdct_NOOS_h)
2525        ENDIF
2526
2527
2528     ENDIF
2529
2530
2531     CALL wrk_dealloc(nb_type , zsumclasses )
2532     
2533     DEALLOCATE(noos_iom_dummy)
2534
2535
2536
2537  END SUBROUTINE dia_dct_wri_NOOS_h
2538
2539  SUBROUTINE dia_dct_wri(kt,ksec,sec)
2540     !!-------------------------------------------------------------
2541     !! Write transport output in numdct
2542     !!
2543     !! Purpose: Write  transports in ascii files
2544     !!
2545     !! Method:
2546     !!        1. Write volume transports in "volume_transport"
2547     !!           Unit: Sv : area * Velocity / 1.e6
2548     !!
2549     !!        2. Write heat transports in "heat_transport"
2550     !!           Unit: Peta W : area * Velocity * T * rhau * Cp / 1.e15
2551     !!
2552     !!        3. Write salt transports in "salt_transport"
2553     !!           Unit: 10^9 g m^3 / s : area * Velocity * S / 1.e6
2554     !!
2555     !!-------------------------------------------------------------
2556     !!arguments
2557     INTEGER, INTENT(IN)          :: kt          ! time-step
2558     TYPE(SECTION), INTENT(INOUT) :: sec         ! section to write   
2559     INTEGER ,INTENT(IN)          :: ksec        ! section number
2560
2561     !!local declarations
2562     INTEGER               :: jclass             ! Dummy loop
2563     CHARACTER(len=2)      :: classe             ! Classname
2564     REAL(wp)              :: zbnd1,zbnd2        ! Class bounds
2565     REAL(wp)              :: zslope             ! section's slope coeff
2566     !
2567     REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace
2568     !!-------------------------------------------------------------
2569     CALL wrk_alloc(nb_type , zsumclasses ) 
2570
2571     zsumclasses(:)=0._wp
2572     zslope = sec%slopeSection       
2573
2574 
2575     DO jclass=1,MAX(1,sec%nb_class-1)
2576
2577        classe   = 'N       '
2578        zbnd1   = 0._wp
2579        zbnd2   = 0._wp
2580        zsumclasses(1:nb_type)=zsumclasses(1:nb_type)+sec%transport(1:nb_type,jclass)
2581
2582   
2583        !insitu density classes transports
2584        IF( ( sec%zsigi(jclass)   .NE. 99._wp ) .AND. &
2585            ( sec%zsigi(jclass+1) .NE. 99._wp )       )THEN
2586           classe = 'DI       '
2587           zbnd1 = sec%zsigi(jclass)
2588           zbnd2 = sec%zsigi(jclass+1)
2589        ENDIF
2590        !potential density classes transports
2591        IF( ( sec%zsigp(jclass)   .NE. 99._wp ) .AND. &
2592            ( sec%zsigp(jclass+1) .NE. 99._wp )       )THEN
2593           classe = 'DP      '
2594           zbnd1 = sec%zsigp(jclass)
2595           zbnd2 = sec%zsigp(jclass+1)
2596        ENDIF
2597        !depth classes transports
2598        IF( ( sec%zlay(jclass)    .NE. 99._wp ) .AND. &
2599            ( sec%zlay(jclass+1)  .NE. 99._wp )       )THEN
2600           classe = 'Z       '
2601           zbnd1 = sec%zlay(jclass)
2602           zbnd2 = sec%zlay(jclass+1)
2603        ENDIF
2604        !salinity classes transports
2605        IF( ( sec%zsal(jclass) .NE. 99._wp    ) .AND. &
2606            ( sec%zsal(jclass+1) .NE. 99._wp  )       )THEN
2607           classe = 'S       '
2608           zbnd1 = sec%zsal(jclass)
2609           zbnd2 = sec%zsal(jclass+1)   
2610        ENDIF
2611        !temperature classes transports
2612        IF( ( sec%ztem(jclass) .NE. 99._wp     ) .AND. &
2613            ( sec%ztem(jclass+1) .NE. 99._wp     )       ) THEN
2614           classe = 'T       '
2615           zbnd1 = sec%ztem(jclass)
2616           zbnd2 = sec%ztem(jclass+1)
2617        ENDIF
2618                 
2619        !write volume transport per class
2620        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, &
2621                              jclass,classe,zbnd1,zbnd2,&
2622                              sec%transport(1,jclass),sec%transport(2,jclass), &
2623                              sec%transport(1,jclass)+sec%transport(2,jclass)
2624
2625        IF( sec%llstrpond )THEN
2626
2627           !write heat transport per class:
2628           WRITE(numdct_temp,119) ndastp,kt,ksec,sec%name,zslope,  &
2629                              jclass,classe,zbnd1,zbnd2,&
2630                              sec%transport(7,jclass)*1000._wp*rcp/1.e15,sec%transport(8,jclass)*1000._wp*rcp/1.e15, &
2631                              ( sec%transport(7,jclass)+sec%transport(8,jclass) )*1000._wp*rcp/1.e15
2632           !write salt transport per class
2633           WRITE(numdct_sal ,119) ndastp,kt,ksec,sec%name,zslope,  &
2634                              jclass,classe,zbnd1,zbnd2,&
2635                              sec%transport(9,jclass)*1000._wp/1.e9,sec%transport(10,jclass)*1000._wp/1.e9,&
2636                              (sec%transport(9,jclass)+sec%transport(10,jclass))*1000._wp/1.e9
2637        ENDIF
2638
2639     ENDDO
2640
2641     zbnd1 = 0._wp
2642     zbnd2 = 0._wp
2643     jclass=0
2644
2645     !write total volume transport
2646     WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, &
2647                           jclass,"total",zbnd1,zbnd2,&
2648                           zsumclasses(1),zsumclasses(2),zsumclasses(1)+zsumclasses(2)
2649
2650     IF( sec%llstrpond )THEN
2651
2652        !write total heat transport
2653        WRITE(numdct_temp,119) ndastp,kt,ksec,sec%name,zslope, &
2654                           jclass,"total",zbnd1,zbnd2,&
2655                           zsumclasses(7)* 1000._wp*rcp/1.e15,zsumclasses(8)* 1000._wp*rcp/1.e15,&
2656                           (zsumclasses(7)+zsumclasses(8) )* 1000._wp*rcp/1.e15
2657        !write total salt transport
2658        WRITE(numdct_sal ,119) ndastp,kt,ksec,sec%name,zslope, &
2659                           jclass,"total",zbnd1,zbnd2,&
2660                           zsumclasses(9)*1000._wp/1.e9,zsumclasses(10)*1000._wp/1.e9,&
2661                           (zsumclasses(9)+zsumclasses(10))*1000._wp/1.e9
2662     ENDIF
2663
2664     
2665     IF ( sec%ll_ice_section) THEN
2666        !write total ice volume transport
2667        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,&
2668                              jclass,"ice_vol",zbnd1,zbnd2,&
2669                              sec%transport(11,1),sec%transport(12,1),&
2670                              sec%transport(11,1)+sec%transport(12,1)
2671        !write total ice surface transport
2672        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,&
2673                              jclass,"ice_surf",zbnd1,zbnd2,&
2674                              sec%transport(13,1),sec%transport(14,1), &
2675                              sec%transport(13,1)+sec%transport(14,1) 
2676     ENDIF
2677                                             
2678118 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3F12.4)
2679119 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6)
2680
2681     CALL wrk_dealloc(nb_type , zsumclasses ) 
2682  END SUBROUTINE dia_dct_wri
2683
2684  FUNCTION interp(ki, kj, kk, cd_point, ptab)
2685  !!----------------------------------------------------------------------
2686  !!
2687  !!   Purpose: compute Temperature/Salinity/density at U-point or V-point
2688  !!   --------
2689  !!
2690  !!   Method:
2691  !!   ------
2692  !!
2693  !!   ====> full step and partial step
2694  !!
2695  !!
2696  !!    |    I          |    I+1           |    Z=Temperature/Salinity/density at U-poinT
2697  !!    |               |                  |
2698  !!  ----------------------------------------  1. Veritcale interpolation: compute zbis
2699  !!    |               |                  |       interpolation between ptab(I,J,K) and ptab(I,J,K+1)
2700  !!    |               |                  |       zbis =
2701  !!    |               |                  |      [ e3w(I+1,J,K)*ptab(I,J,K) + ( e3w(I,J,K) - e3w(I+1,J,K) ) * ptab(I,J,K-1) ]
2702  !!    |               |                  |      /[ e3w(I+1,J,K) + e3w(I,J,K) - e3w(I+1,J,K) ]
2703  !!    |               |                  |
2704  !!    |               |                  |    2. Horizontal interpolation: compute value at U/V point
2705  !!K-1 | ptab(I,J,K-1) |                  |       interpolation between zbis and ptab(I+1,J,K) 
2706  !!    |     .         |                  |
2707  !!    |     .         |                  |       interp = ( 0.5*zet2*zbis + 0.5*zet1*ptab(I+1,J,K) )/(0.5*zet2+0.5*zet1)
2708  !!    |     .         |                  |
2709  !!  ------------------------------------------
2710  !!    |     .         |                  |
2711  !!    |     .         |                  |
2712  !!    |     .         |                  |
2713  !!K   |    zbis.......U...ptab(I+1,J,K)  |
2714  !!    |     .         |                  |
2715  !!    | ptab(I,J,K)   |                  |
2716  !!    |               |------------------|
2717  !!    |               | partials         |
2718  !!    |               |  steps           |
2719  !!  -------------------------------------------
2720  !!    <----zet1------><----zet2--------->
2721  !!
2722  !!
2723  !!   ====>  s-coordinate
2724  !!     
2725  !!    |                |                  |   1. Compute distance between T1 and U points: SQRT( zdep1^2 + (0.5 * zet1 )^2
2726  !!    |                |                  |      Compute distance between T2 and U points: SQRT( zdep2^2 + (0.5 * zet2 )^2
2727  !!    |                | ptab(I+1,J,K)    |
2728  !!    |                |      T2          |   2. Interpolation between  T1 and T2 values at U point
2729  !!    |                |      ^           |   
2730  !!    |                |      | zdep2     |   
2731  !!    |                |      |           |   
2732  !!    |       ^        U      v           |
2733  !!    |       |        |                  |
2734  !!    |       | zdep1  |                  |   
2735  !!    |       v        |                  |
2736  !!    |      T1        |                  |
2737  !!    | ptab(I,J,K)    |                  |
2738  !!    |                |                  |
2739  !!    |                |                  |
2740  !!
2741  !!    <----zet1--------><----zet2--------->
2742  !!
2743  !!----------------------------------------------------------------------
2744  !*arguments
2745  INTEGER, INTENT(IN)                          :: ki, kj, kk   ! coordinate of point
2746  CHARACTER(len=1), INTENT(IN)                 :: cd_point     ! type of point (U, V)
2747  REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ptab         ! variable to compute at (ki, kj, kk )
2748  REAL(wp)                                     :: interp       ! interpolated variable
2749
2750  !*local declations
2751  INTEGER :: ii1, ij1, ii2, ij2                                ! local integer
2752  REAL(wp):: ze3t, zfse3, zwgt1, zwgt2, zbis, zdepu            ! local real
2753  REAL(wp):: zet1, zet2                                        ! weight for interpolation
2754  REAL(wp):: zdep1,zdep2                                       ! differences of depth
2755  !!----------------------------------------------------------------------
2756
2757  IF( cd_point=='U' )THEN
2758     ii1 = ki    ; ij1 = kj 
2759     ii2 = ki+1  ; ij2 = kj 
2760
2761     zet1=e1t(ii1,ij1)
2762     zet2=e1t(ii2,ij2)
2763 
2764
2765  ELSE ! cd_point=='V'
2766     ii1 = ki    ; ij1 = kj 
2767     ii2 = ki    ; ij2 = kj+1 
2768
2769     zet1=e2t(ii1,ij1)
2770     zet2=e2t(ii2,ij2)
2771
2772  ENDIF
2773
2774  IF( ln_sco )THEN   ! s-coordinate case
2775
2776     zdepu = ( fsdept(ii1,ij1,kk) +  fsdept(ii2,ij2,kk) ) /2 
2777     zdep1 = fsdept(ii1,ij1,kk) - zdepu
2778     zdep2 = fsdept(ii2,ij2,kk) - zdepu
2779
2780     !weights
2781     zwgt1 = SQRT( ( 0.5 * zet1 ) * ( 0.5 * zet1 ) + ( zdep1 * zdep1 ) )
2782     zwgt2 = SQRT( ( 0.5 * zet2 ) * ( 0.5 * zet2 ) + ( zdep2 * zdep2 ) )
2783 
2784     ! result
2785     interp = umask(ii1,ij1,kk) * ( zwgt2 *  ptab(ii1,ij1,kk) + zwgt1 *  ptab(ii1,ij1,kk) ) / ( zwgt2 + zwgt1 )   
2786
2787
2788  ELSE       ! full step or partial step case
2789
2790#if defined key_vvl
2791
2792     ze3t  = fse3t_n(ii2,ij2,kk) - fse3t_n(ii1,ij1,kk) 
2793     zwgt1 = ( fse3w_n(ii2,ij2,kk) - fse3w_n(ii1,ij1,kk) ) / fse3w_n(ii2,ij2,kk)
2794     zwgt2 = ( fse3w_n(ii1,ij1,kk) - fse3w_n(ii2,ij2,kk) ) / fse3w_n(ii1,ij1,kk)
2795
2796#else
2797
2798     ze3t  = fse3t(ii2,ij2,kk)   - fse3t(ii1,ij1,kk) 
2799     zwgt1 = ( fse3w(ii2,ij2,kk) - fse3w(ii1,ij1,kk) ) / fse3w(ii2,ij2,kk)
2800     zwgt2 = ( fse3w(ii1,ij1,kk) - fse3w(ii2,ij2,kk) ) / fse3w(ii1,ij1,kk)
2801
2802#endif
2803
2804     IF(kk .NE. 1)THEN
2805
2806        IF( ze3t >= 0. )THEN 
2807           ! zbis
2808           zbis = ptab(ii2,ij2,kk) + zwgt1 * ( ptab(ii2,ij2,kk-1) - ptab(ii2,ij2,kk) ) 
2809           ! result
2810            interp = umask(ii1,ij1,kk) * ( zet2 * ptab(ii1,ij1,kk) + zet1 * zbis )/( zet1 + zet2 )
2811        ELSE
2812           ! zbis
2813           zbis = ptab(ii1,ij1,kk) + zwgt2 * ( ptab(ii1,ij1,kk-1) - ptab(ii1,ij2,kk) )
2814           ! result
2815           interp = umask(ii1,ij1,kk) * ( zet2 * zbis + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 )
2816        ENDIF   
2817
2818     ELSE
2819        interp = umask(ii1,ij1,kk) * (  zet2 * ptab(ii1,ij1,kk) + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 )
2820     ENDIF
2821
2822  ENDIF
2823
2824
2825  END FUNCTION interp
2826
2827#else
2828   !!----------------------------------------------------------------------
2829   !!   Default option :                                       Dummy module
2830   !!----------------------------------------------------------------------
2831   LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .FALSE.    !: diamht flag
2832   PUBLIC 
2833   !! $Id$
2834CONTAINS
2835
2836   SUBROUTINE dia_dct_init          ! Dummy routine
2837      WRITE(*,*) 'dia_dct_init: You should not have seen this print! error?', kt
2838   END SUBROUTINE dia_dct_init
2839
2840   SUBROUTINE dia_dct( kt )         ! Dummy routine
2841      INTEGER, INTENT( in ) :: kt   ! ocean time-step index
2842      WRITE(*,*) 'dia_dct: You should not have seen this print! error?', kt
2843   END SUBROUTINE dia_dct
2844#endif
2845
2846END MODULE diadct
Note: See TracBrowser for help on using the repository browser.