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

Last change on this file since 10492 was 10492, checked in by hadjt, 21 months ago

Updated CO6_shelfclimate_fabm for new NOOS Cross sections

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