New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
diadct.F90 in NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/DIA – NEMO

source: NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/DIA/diadct.F90 @ 15434

Last change on this file since 15434 was 15434, checked in by hadjt, 11 months ago

Region mean and NOOS cross-sections

Use rcp (Specific heat capacity) from nemo physical constant

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