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 @ 15397

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

Region mean and NOOS transport XIOS scaler output
Scalar output, so passing 85 and 3 values per iom call rather than 297*375*85.
DIADCT edited (nd changes commmeted out) to allow compilation on IFORT.
In all NEMO, reading section_ijglobal.diadct was the only compilation error, ifort didn't like ctl_opn, and cray didn't like call ftell.

File size: 140.0 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(5,jsec,jseg,jk) = transports_3d(5,jsec,jseg,jk)  + zTnorm * 0.001 * zsn * 1026._wp
1261!                 transports_3d(5,jsec,jseg,jk) = transports_3d(5,jsec,jseg,jk)  + zTnorm * 0.001 * zsn * zrhop
1262
1263
1264                 transports_3d_inst(2,jsec,jseg,jk) =  zTnorm * zrhoi
1265                 transports_3d_inst(3,jsec,jseg,jk) =  zTnorm * zrhop
1266                 transports_3d_inst(4,jsec,jseg,jk) =  zTnorm * 3850.0 * (ztn) * zrhop ! # 1026._wp !rhop(ji,jj,jk)
1267                 transports_3d_inst(5,jsec,jseg,jk) =  zTnorm * 0.001 * zsn * zrhop
1268
1269
1270                  IF ( ld_update_trans ) THEN
1271                     transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk) + transports_3d_inst(2,jsec,jseg,jk)
1272                     transports_3d(3,jsec,jseg,jk) = transports_3d(3,jsec,jseg,jk) + transports_3d_inst(3,jsec,jseg,jk)
1273                     transports_3d(4,jsec,jseg,jk) = transports_3d(4,jsec,jseg,jk) + transports_3d_inst(4,jsec,jseg,jk)
1274                     transports_3d(5,jsec,jseg,jk) = transports_3d(5,jsec,jseg,jk) + transports_3d_inst(5,jsec,jseg,jk)
1275                  ENDIF
1276
1277
1278
1279
1280              ENDIF
1281   
1282           ENDDO !end of loop on the level
1283
1284#if defined key_lim2 || defined key_lim3
1285
1286           !ICE CASE   
1287           !------------
1288           IF( sec%ll_ice_section )THEN
1289              SELECT CASE (sec%direction(jseg))
1290              CASE(0)
1291                 zumid_ice = 0
1292                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
1293              CASE(1)
1294                 zumid_ice = 0
1295                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
1296              CASE(2)
1297                 zvmid_ice = 0
1298                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
1299              CASE(3)
1300                 zvmid_ice = 0
1301                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
1302              END SELECT
1303   
1304              zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J)
1305   
1306              transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)*   &
1307                                   (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &
1308                                  *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) +  &
1309                                    hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) &
1310                                   +zice_vol_pos
1311              transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)*   &
1312                                    (1.0 -  frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &
1313                                   +zice_surf_pos
1314   
1315           ENDIF !end of ice case
1316#endif
1317 
1318        ENDDO !end of loop on the segment
1319
1320     ENDIF !end of sec%nb_point =0 case
1321     !
1322  END SUBROUTINE transport
1323 
1324  SUBROUTINE transport_h(sec,ld_debug,jsec)
1325     !!-------------------------------------------------------------------------------------------
1326     !!                     ***  ROUTINE hourly_transport  ***
1327     !!
1328     !!              Exactly as routine transport but for data summed at
1329     !!              each time step and averaged each hour
1330     !!
1331     !!  Purpose ::  Compute the transport for each point in a section
1332     !!
1333     !!  Method  ::  Loop over each segment, and each vertical level and add the transport
1334     !!              Be aware :         
1335     !!              One section is a sum of segments
1336     !!              One segment is defined by 2 consecutive points in sec%listPoint
1337     !!              All points of sec%listPoint are positioned on the F-point of the cell
1338     !!
1339     !!              There are two loops:                 
1340     !!              loop on the segment between 2 nodes
1341     !!              loop on the level jk
1342     !!
1343     !! ** Output: Arrays containing the volume,density,salinity,temperature etc
1344     !!            transports for each point in a section, summed over each nn_dct.
1345     !!
1346     !!-------------------------------------------------------------------------------------------
1347     !! * Arguments
1348     TYPE(SECTION),INTENT(INOUT) :: sec
1349     LOGICAL      ,INTENT(IN)    :: ld_debug
1350     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section
1351   
1352     !! * Local variables
1353     INTEGER             :: jk,jseg,jclass,    &!loop on level/segment/classes
1354                            isgnu  , isgnv      !
1355     REAL(wp):: zumid        , zvmid        ,  &!U/V velocity on a cell segment
1356                zumid_ice    , zvmid_ice    ,  &!U/V ice velocity
1357                zTnorm                          !transport of velocity through one cell's sides
1358     REAL(wp):: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point
1359
1360     TYPE(POINT_SECTION) :: k
1361     !!--------------------------------------------------------
1362
1363     !!NIALL IF( ld_debug )WRITE(numout,*)'      Compute transport'
1364
1365     !---------------------------!
1366     !  COMPUTE TRANSPORT        !
1367     !---------------------------!
1368     IF(sec%nb_point .NE. 0)THEN   
1369
1370        !----------------------------------------------------------------------------------------------------
1371        !----------------------------------------------------------------------------------------------------
1372        !----------------------------------------------------------------------------------------------------
1373        !
1374        !
1375        ! ! ! ! JT 1/09/2018 - changing convention. Always direction + is toward left hand of section
1376        !
1377        !    Making sign of the velocities used to calculate the volume transport a function of direction, not slopesection
1378        !    (isgnu, isgnv)
1379        !   
1380        !    They vary for each segment of the section.
1381        !
1382        !----------------------------------------------------------------------------------------------------
1383        !----------------------------------------------------------------------------------------------------
1384        !----------------------------------------------------------------------------------------------------
1385        !Compute sign for velocities:
1386        !
1387        !convention:
1388        !   non horizontal section: direction + is toward left hand of section
1389        !       horizontal section: direction + is toward north of section
1390        !
1391        !
1392        !       slopeSection < 0     slopeSection > 0       slopeSection=inf            slopeSection=0
1393        !       ----------------      -----------------     ---------------             --------------
1394        !
1395        !   isgnv=1         direction +     
1396        !  ______         _____             ______                                                   
1397        !        |           //|            |                  |                         direction +   
1398        !        | isgnu=1  // |            |isgnu=1           |isgnu=1                     /|\
1399        !        |_______  //         ______|    \\            | ---\                        |
1400        !               |             | isgnv=-1  \\ |         | ---/ direction +       ____________
1401        !               |             |          __\\|         |                   
1402        !               |             |     direction +        |                      isgnv=1                                 
1403        !                                                     
1404        !----------------------------------------------------------------------------------------------------
1405
1406        IF( ld_debug )write(numout,*)"isgnu isgnv ",isgnu,isgnv
1407
1408        !--------------------------------------!
1409        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  !
1410        !--------------------------------------!
1411        DO jseg=1,MAX(sec%nb_point-1,0)
1412           
1413           
1414           !Compute sign for velocities:
1415           
1416           !isgnu =  1
1417           !isgnv =  1
1418           !
1419           ! changing sign of u and v is dependent on the direction of the section.
1420           !isgnu =  1
1421           !isgnv =  1
1422           !SELECT CASE( sec%direction(jseg) )
1423           !CASE(0)  ;   isgnv = -1
1424           !CASE(3)  ;   isgnu = -1
1425           !END SELECT
1426           
1427           
1428           SELECT CASE( sec%direction(jseg) )
1429           CASE(0) 
1430              isgnu =  1
1431              isgnv = -1
1432           CASE(1)
1433              isgnu =  1
1434              isgnv =  1
1435           CASE(2) 
1436              isgnu =  1
1437              isgnv =  1
1438           CASE(3) 
1439              isgnu = -1
1440              isgnv =  1
1441           END SELECT
1442           
1443           !-------------------------------------------------------------------------------------------
1444           ! Select the appropriate coordinate for computing the velocity of the segment
1445           ! Corrected by JT 01/09/2018 (#)
1446           !
1447           !                      CASE(0)                                    Case (2)
1448           !                      -------                                    --------
1449           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)     
1450           !      F(i,j)---------#V(i,j)-------F(i+1,j)                                 |
1451           !                     -------->                                              |
1452           !                                                                        |   |
1453           !                                                                        |   |
1454           !                      Case (3)                                          | U(i,j)
1455           !                      --------                                          |   |
1456           !                                                                        V   |
1457           !  listPoint(jseg+1) F(i,j+1)                                                |
1458           !                        |                                                   |
1459           !                        |                                                   |
1460           !                        |                                 listPoint(jseg+1) F(i,j-1)
1461           !                   ^    |                                           
1462           !                   |    |                                           
1463           !                   | U(i,j+1)                                           
1464           !                   |    |                                       Case(1)     
1465           !                   |    |                                       ------     
1466           !                        |                                           
1467           !                        |                 listPoint(jseg+1)             listPoint(jseg)                           
1468           !                        |                 F(i-1,j)----------#V(i-1,j) ------#f(i,j)                           
1469           ! listPoint(jseg)     F(i,j)                                 <-------
1470           !
1471           !-------------------------------------------------------------------------------------------
1472
1473           SELECT CASE( sec%direction(jseg) )
1474           CASE(0)  ;   k = sec%listPoint(jseg)
1475           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
1476           CASE(2)  ;   k = sec%listPoint(jseg)
1477           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
1478           END SELECT
1479
1480           !---------------------------|
1481           !     LOOP ON THE LEVEL     |
1482           !---------------------------|
1483           !Sum of the transport on the vertical
1484           DO jk=1,mbkt(k%I,k%J)   ! mbathy(k%I,k%J)
1485
1486              ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
1487              SELECT CASE( sec%direction(jseg) )
1488              CASE(0,1)
1489                 ztn   = interp(k%I,k%J,jk,'V',0)
1490                 zsn   = interp(k%I,k%J,jk,'V',1)
1491                 zrhop = interp(k%I,k%J,jk,'V',2)
1492                 zrhoi = interp(k%I,k%J,jk,'V',3)
1493                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1)
1494              CASE(2,3)
1495                 ztn   = interp(k%I,k%J,jk,'U',0)
1496                 zsn   = interp(k%I,k%J,jk,'U',1)
1497                 zrhop = interp(k%I,k%J,jk,'U',2)
1498                 zrhoi = interp(k%I,k%J,jk,'U',3)
1499                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1) 
1500              END SELECT
1501
1502              zfsdep = gdept_n(k%I,k%J,jk)
1503 
1504              !compute velocity with the correct direction
1505              SELECT CASE( sec%direction(jseg) )
1506              CASE(0,1) 
1507                 zumid=0.
1508                 zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk)
1509              CASE(2,3)
1510                 zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk)
1511                 zvmid=0.
1512              END SELECT
1513
1514              !zTnorm=transport through one cell;
1515              !velocity* cell's length * cell's thickness
1516              !zTnorm=zumid*e2u(k%I,k%J)*  fse3u(k%I,k%J,jk)+     &
1517              !       zvmid*e1v(k%I,k%J)*  fse3v(k%I,k%J,jk)
1518              zTnorm=zumid*e2u(k%I,k%J)*  e3u_n(k%I,k%J,jk)+     &
1519                     zvmid*e1v(k%I,k%J)*  e3v_n(k%I,k%J,jk)
1520
1521#if ! defined key_vvl
1522              !add transport due to free surface
1523              IF( jk==1 )THEN
1524                 zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + &
1525                                   zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk)
1526              ENDIF
1527#endif
1528              !COMPUTE TRANSPORT
1529
1530              transports_3d_h(1,jsec,jseg,jk) = transports_3d_h(1,jsec,jseg,jk) + zTnorm
1531 
1532              IF ( sec%llstrpond ) THEN
1533                 transports_3d_h(2,jsec,jseg,jk) = transports_3d_h(2,jsec,jseg,jk)  + zTnorm * zrhoi
1534                 transports_3d_h(3,jsec,jseg,jk) = transports_3d_h(3,jsec,jseg,jk)  + zTnorm * zrhop
1535                 !transports_3d_h(4,jsec,jseg,jk) = transports_3d_h(4,jsec,jseg,jk)  + zTnorm * 4.e+3_wp * (ztn+273.15) * 1026._wp
1536                 transports_3d_h(4,jsec,jseg,jk) = transports_3d_h(4,jsec,jseg,jk)  + zTnorm * 3850.0 * (ztn) * zrhop
1537                 !transports_3d_h(5,jsec,jseg,jk) = transports_3d_h(5,jsec,jseg,jk)  + zTnorm * 0.001 * zsn * 1026._wp
1538                 transports_3d_h(5,jsec,jseg,jk) = transports_3d_h(5,jsec,jseg,jk)  + zTnorm * 0.001 * zsn * zrhop
1539              ENDIF
1540
1541           ENDDO !end of loop on the level
1542
1543#if defined key_lim2 || defined key_lim3
1544
1545           !ICE CASE   
1546           !------------
1547           IF( sec%ll_ice_section )THEN
1548              SELECT CASE (sec%direction(jseg))
1549              CASE(0)
1550                 zumid_ice = 0
1551                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
1552              CASE(1)
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(2)
1556                 zvmid_ice = 0
1557                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
1558              CASE(3)
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              END SELECT
1562   
1563              zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J)
1564   
1565              transports_2d_h(1,jsec,jseg) = transports_2d_h(1,jsec,jseg) + (zTnorm)*   &
1566                                   (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &
1567                                  *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) +  &
1568                                    hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) &
1569                                   +zice_vol_pos
1570              transports_2d_h(2,jsec,jseg) = transports_2d_h(2,jsec,jseg) + (zTnorm)*   &
1571                                    (1.0 -  frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &
1572                                   +zice_surf_pos
1573   
1574           ENDIF !end of ice case
1575#endif
1576 
1577        ENDDO !end of loop on the segment
1578
1579     ENDIF   !end of sec%nb_point =0 case
1580     !
1581  END SUBROUTINE transport_h
1582 
1583  SUBROUTINE dia_dct_sum(sec,jsec) 
1584     !!-------------------------------------------------------------
1585     !! Purpose: Average the transport over nn_dctwri time steps 
1586     !! and sum over the density/salinity/temperature/depth classes
1587     !!
1588     !! Method:   Sum over relevant grid cells to obtain values 
1589     !!           for each class
1590     !!              There are several loops:                 
1591     !!              loop on the segment between 2 nodes
1592     !!              loop on the level jk
1593     !!              loop on the density/temperature/salinity/level classes
1594     !!              test on the density/temperature/salinity/level
1595     !!
1596     !!  Note:    Transport through a given section is equal to the sum of transports
1597     !!           computed on each proc.
1598     !!           On each proc,transport is equal to the sum of transport computed through
1599     !!           segments linking each point of sec%listPoint  with the next one.   
1600     !!
1601     !!-------------------------------------------------------------
1602     !! * arguments
1603     TYPE(SECTION),INTENT(INOUT) :: sec 
1604     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section
1605 
1606     TYPE(POINT_SECTION) :: k 
1607     INTEGER  :: jk,jseg,jclass                        ! dummy variables for looping on level/segment/classes 
1608     REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point
1609     !!-------------------------------------------------------------
1610
1611
1612     !! Sum the relevant segments to obtain values for each class
1613     IF(sec%nb_point .NE. 0)THEN   
1614
1615        !--------------------------------------!
1616        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  !
1617        !--------------------------------------!
1618        DO jseg=1,MAX(sec%nb_point-1,0)
1619           
1620           !-------------------------------------------------------------------------------------------
1621           ! Select the appropriate coordinate for computing the velocity of the segment
1622           !
1623           !                      CASE(0)                                    Case (2)
1624           !                      -------                                    --------
1625           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)     
1626           !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               |
1627           !                                                                            |
1628           !                                                                            |
1629           !                                                                            |
1630           !                      Case (3)                                            U(i,j)
1631           !                      --------                                              |
1632           !                                                                            |
1633           !  listPoint(jseg+1) F(i,j+1)                                                |
1634           !                        |                                                   |
1635           !                        |                                                   |
1636           !                        |                                 listPoint(jseg+1) F(i,j-1)
1637           !                        |                                           
1638           !                        |                                           
1639           !                     U(i,j+1)                                           
1640           !                        |                                       Case(1)     
1641           !                        |                                       ------     
1642           !                        |                                           
1643           !                        |                 listPoint(jseg+1)             listPoint(jseg)                           
1644           !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                           
1645           ! listPoint(jseg)     F(i,j)
1646           !
1647           !-------------------------------------------------------------------------------------------
1648
1649           SELECT CASE( sec%direction(jseg) )
1650           CASE(0)  ;   k = sec%listPoint(jseg)
1651           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
1652           CASE(2)  ;   k = sec%listPoint(jseg)
1653           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
1654           END SELECT
1655
1656           !---------------------------|
1657           !     LOOP ON THE LEVEL     |
1658           !---------------------------|
1659           !Sum of the transport on the vertical
1660           DO jk=1,mbkt(k%I,k%J)   ! mbathy(k%I,k%J)
1661
1662              ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
1663              SELECT CASE( sec%direction(jseg) )
1664              CASE(0,1)
1665                 ztn   = interp(k%I,k%J,jk,'V',0)
1666                 zsn   = interp(k%I,k%J,jk,'V',1)
1667                 zrhop = interp(k%I,k%J,jk,'V',2)
1668                 zrhoi = interp(k%I,k%J,jk,'V',3)
1669                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1)
1670              CASE(2,3)
1671                 ztn   = interp(k%I,k%J,jk,'U',0)
1672                 zsn   = interp(k%I,k%J,jk,'U',1)
1673                 zrhop = interp(k%I,k%J,jk,'U',2)
1674                 zrhoi = interp(k%I,k%J,jk,'U',3)
1675                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1) 
1676              END SELECT
1677
1678              zfsdep = gdept_n(k%I,k%J,jk) 
1679 
1680              !-------------------------------
1681              !  LOOP ON THE DENSITY CLASSES |
1682              !-------------------------------
1683              !The computation is made for each density/temperature/salinity/depth class
1684              DO jclass=1,MAX(1,sec%nb_class-1)
1685 
1686                 !----------------------------------------------!
1687                 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL!
1688                 !----------------------------------------------!
1689 
1690                 IF ( (                                                    &
1691                    ((( zrhop .GE. (sec%zsigp(jclass)+1000.  )) .AND.      &
1692                    (   zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR.     &
1693                    ( sec%zsigp(jclass) .EQ. 99.)) .AND.                   &
1694
1695                    ((( zrhoi .GE. (sec%zsigi(jclass) + 1000.  )) .AND.    &
1696                    (   zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR.     &
1697                    ( sec%zsigi(jclass) .EQ. 99.)) .AND.                   &
1698
1699                    ((( zsn .GT. sec%zsal(jclass)) .AND.                   &
1700                    (   zsn .LE. sec%zsal(jclass+1))) .OR.                 &
1701                    ( sec%zsal(jclass) .EQ. 99.)) .AND.                    &
1702
1703                    ((( ztn .GE. sec%ztem(jclass)) .AND.                   &
1704                    (   ztn .LE. sec%ztem(jclass+1))) .OR.                 &
1705                    ( sec%ztem(jclass) .EQ.99.)) .AND.                     &
1706
1707                    ((( zfsdep .GE. sec%zlay(jclass)) .AND.                &
1708                    (   zfsdep .LE. sec%zlay(jclass+1))) .OR.              &
1709                    ( sec%zlay(jclass) .EQ. 99. ))                         &
1710                                                                   ))   THEN
1711
1712                    !SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS
1713                    !----------------------------------------------------------------------------
1714                    IF (transports_3d(1,jsec,jseg,jk) .GE. 0.0) THEN
1715                       sec%transport(1,jclass) = sec%transport(1,jclass)+transports_3d(1,jsec,jseg,jk)
1716                    ELSE
1717                       sec%transport(2,jclass) = sec%transport(2,jclass)+transports_3d(1,jsec,jseg,jk)
1718                    ENDIF
1719                    IF( sec%llstrpond )THEN
1720
1721                       IF( transports_3d(1,jsec,jseg,jk) .NE. 0._wp ) THEN
1722
1723                          IF (transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) .GE. 0.0 ) THEN
1724                             sec%transport(3,jclass) = sec%transport(3,jclass)+transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk)
1725                          ELSE
1726                             sec%transport(4,jclass) = sec%transport(4,jclass)+transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk)
1727                          ENDIF
1728
1729                          IF ( transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) .GE. 0.0 ) THEN
1730                             sec%transport(5,jclass) = sec%transport(5,jclass)+transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk)
1731                          ELSE
1732                             sec%transport(6,jclass) = sec%transport(6,jclass)+transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk)
1733                          ENDIF
1734
1735                       ENDIF
1736
1737                       IF ( transports_3d(4,jsec,jseg,jk) .GE. 0.0 ) THEN
1738                          sec%transport(7,jclass) = sec%transport(7,jclass)+transports_3d(4,jsec,jseg,jk)
1739                       ELSE
1740                          sec%transport(8,jclass) = sec%transport(8,jclass)+transports_3d(4,jsec,jseg,jk)
1741                       ENDIF
1742
1743                       IF ( transports_3d(5,jsec,jseg,jk) .GE. 0.0 ) THEN
1744                          sec%transport( 9,jclass) = sec%transport( 9,jclass)+transports_3d(5,jsec,jseg,jk)
1745                       ELSE
1746                          sec%transport(10,jclass) = sec%transport(10,jclass)+transports_3d(5,jsec,jseg,jk)
1747                       ENDIF
1748 
1749                    ELSE
1750                       sec%transport( 3,jclass) = 0._wp 
1751                       sec%transport( 4,jclass) = 0._wp 
1752                       sec%transport( 5,jclass) = 0._wp 
1753                       sec%transport( 6,jclass) = 0._wp 
1754                       sec%transport( 7,jclass) = 0._wp
1755                       sec%transport( 8,jclass) = 0._wp
1756                       sec%transport( 9,jclass) = 0._wp
1757                       sec%transport(10,jclass) = 0._wp
1758                    ENDIF
1759 
1760                 ENDIF ! end of test if point is in class
1761   
1762              ENDDO ! end of loop on the classes
1763 
1764           ENDDO ! loop over jk
1765 
1766#if defined key_lim2 || defined key_lim3 
1767 
1768           !ICE CASE     
1769           IF( sec%ll_ice_section )THEN
1770 
1771              IF ( transports_2d(1,jsec,jseg) .GE. 0.0 ) THEN
1772                 sec%transport(11,1) = sec%transport(11,1)+transports_2d(1,jsec,jseg)
1773              ELSE
1774                 sec%transport(12,1) = sec%transport(12,1)+transports_2d(1,jsec,jseg)
1775              ENDIF
1776 
1777              IF ( transports_2d(3,jsec,jseg) .GE. 0.0 ) THEN
1778                 sec%transport(13,1) = sec%transport(13,1)+transports_2d(2,jsec,jseg)
1779              ELSE
1780                 sec%transport(14,1) = sec%transport(14,1)+transports_2d(2,jsec,jseg)
1781              ENDIF
1782 
1783           ENDIF !end of ice case
1784#endif
1785 
1786        ENDDO !end of loop on the segment
1787
1788     ELSE  !if sec%nb_point =0
1789        sec%transport(1:2,:)=0.
1790        IF (sec%llstrpond) sec%transport(3:10,:)=0.
1791        IF (sec%ll_ice_section) sec%transport( 11:14,:)=0.
1792     ENDIF !end of sec%nb_point =0 case
1793
1794  END SUBROUTINE dia_dct_sum
1795 
1796  SUBROUTINE dia_dct_sum_h(sec,jsec)
1797     !!-------------------------------------------------------------
1798     !! Exactly as dia_dct_sum but for hourly files containing data summed at each time step
1799     !!
1800     !! Purpose: Average the transport over nn_dctwri time steps
1801     !! and sum over the density/salinity/temperature/depth classes
1802     !!
1803     !! Method:
1804     !!           Sum over relevant grid cells to obtain values
1805     !!           for each
1806     !!              There are several loops:                 
1807     !!              loop on the segment between 2 nodes
1808     !!              loop on the level jk
1809     !!              loop on the density/temperature/salinity/level classes
1810     !!              test on the density/temperature/salinity/level
1811     !!
1812     !!  ** Method  :Transport through a given section is equal to the sum of transports
1813     !!              computed on each proc.
1814     !!              On each proc,transport is equal to the sum of transport computed through
1815     !!              segments linking each point of sec%listPoint  with the next one.   
1816     !!
1817     !!-------------------------------------------------------------
1818     !! * arguments
1819     TYPE(SECTION),INTENT(INOUT) :: sec
1820     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section
1821
1822     TYPE(POINT_SECTION) :: k
1823     INTEGER  :: jk,jseg,jclass                        !loop on level/segment/classes
1824     REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point
1825     !!-------------------------------------------------------------
1826
1827     !! Sum the relevant segments to obtain values for each class
1828     IF(sec%nb_point .NE. 0)THEN   
1829
1830        !--------------------------------------!
1831        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  !
1832        !--------------------------------------!
1833        DO jseg=1,MAX(sec%nb_point-1,0)
1834           
1835           !-------------------------------------------------------------------------------------------
1836           ! Select the appropriate coordinate for computing the velocity of the segment
1837           !
1838           !                      CASE(0)                                    Case (2)
1839           !                      -------                                    --------
1840           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)     
1841           !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               |
1842           !                                                                            |
1843           !                                                                            |
1844           !                                                                            |
1845           !                      Case (3)                                            U(i,j)
1846           !                      --------                                              |
1847           !                                                                            |
1848           !  listPoint(jseg+1) F(i,j+1)                                                |
1849           !                        |                                                   |
1850           !                        |                                                   |
1851           !                        |                                 listPoint(jseg+1) F(i,j-1)
1852           !                        |                                           
1853           !                        |                                           
1854           !                     U(i,j+1)                                           
1855           !                        |                                       Case(1)     
1856           !                        |                                       ------     
1857           !                        |                                           
1858           !                        |                 listPoint(jseg+1)             listPoint(jseg)                           
1859           !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                           
1860           ! listPoint(jseg)     F(i,j)
1861           !
1862           !-------------------------------------------------------------------------------------------
1863
1864           SELECT CASE( sec%direction(jseg) )
1865           CASE(0)  ;   k = sec%listPoint(jseg)
1866           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
1867           CASE(2)  ;   k = sec%listPoint(jseg)
1868           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
1869           END SELECT
1870
1871           !---------------------------|
1872           !     LOOP ON THE LEVEL     |
1873           !---------------------------|
1874           !Sum of the transport on the vertical
1875           DO jk=1,mbkt(k%I,k%J)   ! mbathy(k%I,k%J)
1876
1877              ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
1878              SELECT CASE( sec%direction(jseg) )
1879              CASE(0,1)
1880                 ztn   = interp(k%I,k%J,jk,'V',0)
1881                 zsn   = interp(k%I,k%J,jk,'V',1)
1882                 zrhop = interp(k%I,k%J,jk,'V',2)
1883                 zrhoi = interp(k%I,k%J,jk,'V',3)
1884                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1)
1885              CASE(2,3)
1886                 ztn   = interp(k%I,k%J,jk,'U',0)
1887                 zsn   = interp(k%I,k%J,jk,'U',1)
1888                 zrhop = interp(k%I,k%J,jk,'U',2)
1889                 zrhoi = interp(k%I,k%J,jk,'U',3)
1890                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1) 
1891              END SELECT
1892
1893              zfsdep = gdept_n(k%I,k%J,jk)
1894 
1895              !-------------------------------
1896              !  LOOP ON THE DENSITY CLASSES |
1897              !-------------------------------
1898              !The computation is made for each density/heat/salt/... class
1899              DO jclass=1,MAX(1,sec%nb_class-1)
1900
1901                 !----------------------------------------------!
1902                 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL!
1903                 !----------------------------------------------!
1904 
1905                 IF ( (                                                    &
1906                    ((( zrhop .GE. (sec%zsigp(jclass)+1000.  )) .AND.      &
1907                    (   zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR.     &
1908                    ( sec%zsigp(jclass) .EQ. 99.)) .AND.                   &
1909
1910                    ((( zrhoi .GE. (sec%zsigi(jclass) + 1000.  )) .AND.    &
1911                    (   zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR.     &
1912                    ( sec%zsigi(jclass) .EQ. 99.)) .AND.                   &
1913
1914                    ((( zsn .GT. sec%zsal(jclass)) .AND.                   &
1915                    (   zsn .LE. sec%zsal(jclass+1))) .OR.                 &
1916                    ( sec%zsal(jclass) .EQ. 99.)) .AND.                    &
1917
1918                    ((( ztn .GE. sec%ztem(jclass)) .AND.                   &
1919                    (   ztn .LE. sec%ztem(jclass+1))) .OR.                 &
1920                    ( sec%ztem(jclass) .EQ.99.)) .AND.                     &
1921
1922                    ((( zfsdep .GE. sec%zlay(jclass)) .AND.                &
1923                    (   zfsdep .LE. sec%zlay(jclass+1))) .OR.              &
1924                    ( sec%zlay(jclass) .EQ. 99. ))                         &
1925                                                                   ))   THEN
1926
1927                    !SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS
1928                    !----------------------------------------------------------------------------
1929                    IF (transports_3d_h(1,jsec,jseg,jk) .GE. 0.0) THEN
1930                       sec%transport_h(1,jclass) = sec%transport_h(1,jclass)+transports_3d_h(1,jsec,jseg,jk)
1931                    ELSE
1932                       sec%transport_h(2,jclass) = sec%transport_h(2,jclass)+transports_3d_h(1,jsec,jseg,jk)
1933                    ENDIF
1934                    IF( sec%llstrpond )THEN
1935
1936                       IF( transports_3d_h(1,jsec,jseg,jk) .NE. 0._wp ) THEN
1937
1938                          IF (transports_3d_h(2,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) .GE. 0.0 ) THEN
1939                             sec%transport_h(3,jclass) = sec%transport_h(3,jclass)+transports_3d_h(2,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk)
1940                          ELSE
1941                             sec%transport_h(4,jclass) = sec%transport_h(4,jclass)+transports_3d_h(2,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk)
1942                          ENDIF
1943
1944                          IF ( transports_3d_h(3,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) .GE. 0.0 ) THEN
1945                             sec%transport_h(5,jclass) = sec%transport_h(5,jclass)+transports_3d_h(3,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk)
1946                          ELSE
1947                             sec%transport_h(6,jclass) = sec%transport_h(6,jclass)+transports_3d_h(3,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk)
1948                          ENDIF
1949
1950                       ENDIF
1951
1952                       IF ( transports_3d_h(4,jsec,jseg,jk) .GE. 0.0 ) THEN
1953                          sec%transport_h(7,jclass) = sec%transport_h(7,jclass)+transports_3d_h(4,jsec,jseg,jk)
1954                       ELSE
1955                          sec%transport_h(8,jclass) = sec%transport_h(8,jclass)+transports_3d_h(4,jsec,jseg,jk)
1956                       ENDIF
1957
1958                       IF ( transports_3d_h(5,jsec,jseg,jk) .GE. 0.0 ) THEN
1959                          sec%transport_h( 9,jclass) = sec%transport_h( 9,jclass)+transports_3d_h(5,jsec,jseg,jk)
1960                       ELSE
1961                          sec%transport_h(10,jclass) = sec%transport_h(10,jclass)+transports_3d_h(5,jsec,jseg,jk)
1962                       ENDIF
1963
1964                    ELSE
1965                       sec%transport_h( 3,jclass) = 0._wp
1966                       sec%transport_h( 4,jclass) = 0._wp
1967                       sec%transport_h( 5,jclass) = 0._wp
1968                       sec%transport_h( 6,jclass) = 0._wp
1969                       sec%transport_h( 7,jclass) = 0._wp
1970                       sec%transport_h( 8,jclass) = 0._wp
1971                       sec%transport_h( 9,jclass) = 0._wp
1972                       sec%transport_h(10,jclass) = 0._wp
1973                    ENDIF
1974
1975                 ENDIF ! end of test if point is in class
1976   
1977              ENDDO ! end of loop on the classes
1978
1979           ENDDO ! loop over jk
1980
1981#if defined key_lim2 || defined key_lim3
1982
1983           !ICE CASE   
1984           IF( sec%ll_ice_section )THEN
1985
1986              IF ( transports_2d_h(1,jsec,jseg) .GE. 0.0 ) THEN
1987                 sec%transport_h(11,1) = sec%transport_h(11,1)+transports_2d_h(1,jsec,jseg)
1988              ELSE
1989                 sec%transport_h(12,1) = sec%transport_h(12,1)+transports_2d_h(1,jsec,jseg)
1990              ENDIF
1991
1992              IF ( transports_2d_h(3,jsec,jseg) .GE. 0.0 ) THEN
1993                 sec%transport_h(13,1) = sec%transport_h(13,1)+transports_2d_h(2,jsec,jseg)
1994              ELSE
1995                 sec%transport_h(14,1) = sec%transport_h(14,1)+transports_2d_h(2,jsec,jseg)
1996              ENDIF
1997
1998           ENDIF !end of ice case
1999#endif
2000 
2001        ENDDO !end of loop on the segment
2002
2003     ELSE  !if sec%nb_point =0
2004        sec%transport_h(1:2,:)=0.
2005        IF (sec%llstrpond) sec%transport_h(3:10,:)=0.
2006        IF (sec%ll_ice_section) sec%transport_h( 11:14,:)=0.
2007     ENDIF !end of sec%nb_point =0 case
2008
2009  END SUBROUTINE dia_dct_sum_h
2010 
2011  SUBROUTINE dia_dct_wri_NOOS(kt,ksec,sec)
2012     !!-------------------------------------------------------------
2013     !! Write transport output in numdct using NOOS formatting
2014     !!
2015     !! Purpose: Write  transports in ascii files
2016     !!
2017     !! Method:
2018     !!        1. Write volume transports in "volume_transport"
2019     !!           Unit: Sv : area * Velocity / 1.e6
2020     !!
2021     !!        2. Write heat transports in "heat_transport"
2022     !!           Unit: Peta W : area * Velocity * T * rhau * Cp / 1.e15
2023     !!
2024     !!        3. Write salt transports in "salt_transport"
2025     !!           Unit: 10^9 g m^3 / s : area * Velocity * S / 1.e6
2026     !!
2027     !!-------------------------------------------------------------
2028     !!arguments
2029     INTEGER, INTENT(IN)          :: kt          ! time-step
2030     TYPE(SECTION), INTENT(INOUT) :: sec         ! section to write   
2031     INTEGER ,INTENT(IN)          :: ksec        ! section number
2032
2033     !!local declarations
2034     INTEGER               :: jclass,ji             ! Dummy loop
2035     CHARACTER(len=2)      :: classe             ! Classname
2036     REAL(wp)              :: zbnd1,zbnd2        ! Class bounds
2037     REAL(wp)              :: zslope             ! section's slope coeff
2038     !
2039     REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace
2040     CHARACTER(len=3)      :: noos_sect_name             ! Classname
2041     CHARACTER(len=25)      :: noos_var_sect_name
2042     REAL(wp), ALLOCATABLE, DIMENSION(:) ::   noos_iom_dummy
2043     INTEGER               :: IERR
2044     
2045     REAL(wp), DIMENSION(3) :: tmp_iom_output
2046     REAL(wp)               :: max_iom_val
2047     LOGICAL       ::   verbose     
2048     verbose = ln_dct_verbose! .false.
2049     
2050     !!-------------------------------------------------------------
2051     
2052     
2053     
2054     IF( lwp .AND. verbose ) THEN
2055        WRITE(numout,*) " "
2056        WRITE(numout,*) "dia_dct_wri_NOOS: write transports through sections at timestep: ", kt
2057        WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
2058     ENDIF   
2059       
2060     !CALL wrk_alloc(nb_type_class , zsumclasses ) 
2061     ALLOCATE( zsumclasses(nb_type_class),  STAT= ierr )
2062     IF( ierr /= 0 )   CALL ctl_stop( 'dia_dct_wri_NOOS: failed to allocate zsumclasses array' )
2063
2064     zsumclasses(:)=0._wp
2065     zslope = sec%slopeSection       
2066     
2067     IF( lwp ) THEN
2068         IF  ( ln_dct_ascii ) THEN
2069             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
2070         ELSE
2071             WRITE(numdct_NOOS) nyear,nmonth,nday,ksec-1,sec%nb_class-1,sec%name
2072         ENDIF 
2073     ENDIF
2074   
2075     ! Sum all classes together, to give one values per type (pos tran, neg vol trans etc...).
2076     DO jclass=1,MAX(1,sec%nb_class-1)
2077        zsumclasses(1:nb_type_class)=zsumclasses(1:nb_type_class)+sec%transport(1:nb_type_class,jclass)
2078     ENDDO
2079 
2080     classe   = 'total   '
2081     zbnd1   = 0._wp
2082     zbnd2   = 0._wp
2083     
2084     
2085     
2086     write (noos_sect_name, "(I0.3)")  ksec
2087     
2088     IF ( nn_dct_iom_cont  .eq. 1) THEN
2089         max_iom_val = 1.e10
2090         ALLOCATE( noos_iom_dummy(3),  STAT= ierr )
2091            IF( ierr /= 0 )   CALL ctl_stop( 'dia_dct_wri_NOOS: failed to allocate noos_iom_dummy array' )
2092     ENDIF
2093     
2094!     JT   I think changing the sign on the output based on the zslope value is redunant.
2095!
2096!     IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN
2097!         
2098!         IF( lwp ) THEN
2099!             WRITE(numdct_NOOS,'(9e12.4E2)') -(zsumclasses( 1)+zsumclasses( 2)), -zsumclasses( 2),-zsumclasses( 1),   &
2100!                                            -(zsumclasses( 7)+zsumclasses( 8)), -zsumclasses( 8),-zsumclasses( 7),   &
2101!                                            -(zsumclasses( 9)+zsumclasses(10)), -zsumclasses(10),-zsumclasses( 9)
2102!             CALL FLUSH(numdct_NOOS)
2103!         endif
2104
2105
2106
2107
2108
2109
2110
2111     
2112     
2113     IF( lwp ) THEN
2114        IF  ( ln_dct_ascii ) THEN
2115             !WRITE(numdct_NOOS,'(9e12.4E2)')   zsumclasses( 1)+zsumclasses( 2) ,  zsumclasses( 1), zsumclasses( 2),   &
2116             WRITE(numdct_NOOS,'(3F18.3,6e16.8E2)')   zsumclasses( 1)+zsumclasses( 2) ,  zsumclasses( 1), zsumclasses( 2),   &
2117                                             zsumclasses( 7)+zsumclasses( 8) ,  zsumclasses( 7), zsumclasses( 8),   &
2118                                             zsumclasses( 9)+zsumclasses(10) ,  zsumclasses( 9), zsumclasses(10)
2119             CALL FLUSH(numdct_NOOS)
2120        ELSE
2121             WRITE(numdct_NOOS)   zsumclasses( 1)+zsumclasses( 2) ,  zsumclasses( 1), zsumclasses( 2),   &
2122                                  zsumclasses( 7)+zsumclasses( 8) ,  zsumclasses( 7), zsumclasses( 8),   &
2123                                  zsumclasses( 9)+zsumclasses(10) ,  zsumclasses( 9), zsumclasses(10)
2124             CALL FLUSH(numdct_NOOS) 
2125         ENDIF
2126     ENDIF
2127!         
2128    IF ( nn_dct_iom_cont .EQ. 1) THEN
2129         noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_trans'
2130         IF (iom_use(noos_var_sect_name)) THEN
2131             noos_iom_dummy(:) = 0.
2132             tmp_iom_output(:) = 0.
2133             
2134             tmp_iom_output(1) = (zsumclasses( 1)+zsumclasses( 2))
2135             tmp_iom_output(2) = zsumclasses( 1)
2136             tmp_iom_output(3) = zsumclasses( 2)
2137             
2138             ! Convert to Sv
2139             tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
2140             tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
2141             tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
2142             
2143             ! limit maximum and minimum values in iom_put
2144             if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
2145             if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
2146             if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
2147             if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
2148             if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
2149             if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
2150             
2151             ! Set NaN's to Zero         
2152             if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
2153             if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
2154             if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
2155             
2156             noos_iom_dummy(1) = tmp_iom_output(1)
2157             noos_iom_dummy(2) = tmp_iom_output(2)
2158             noos_iom_dummy(3) = tmp_iom_output(3)
2159             
2160             !noos_iom_dummy(1) = (zsumclasses( 1)+zsumclasses( 2))
2161             !noos_iom_dummy(2) = zsumclasses( 1)
2162             !noos_iom_dummy(3) = zsumclasses( 2)
2163             
2164             
2165             
2166             if ( lwp .AND. verbose ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
2167             CALL iom_put( noos_var_sect_name,  noos_iom_dummy(:) )
2168         ENDIF
2169
2170         noos_var_sect_name =  "noos_" // trim(noos_sect_name) // '_heat'
2171         IF (iom_use(noos_var_sect_name)) THEN
2172             noos_iom_dummy(:) = 0.
2173             tmp_iom_output(:) = 0.
2174             
2175             tmp_iom_output(1) = (zsumclasses( 7)+zsumclasses( 8))
2176             tmp_iom_output(2) = zsumclasses( 7)
2177             tmp_iom_output(3) = zsumclasses( 8)
2178             
2179             ! Convert to TJ/s
2180             tmp_iom_output(1) = tmp_iom_output(1)*1.E-12
2181             tmp_iom_output(2) = tmp_iom_output(2)*1.E-12
2182             tmp_iom_output(3) = tmp_iom_output(3)*1.E-12
2183             
2184             ! limit maximum and minimum values in iom_put
2185             if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
2186             if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
2187             if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
2188             if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
2189             if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
2190             if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
2191             
2192             ! Set NaN's to Zero         
2193             if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
2194             if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
2195             if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
2196             
2197             noos_iom_dummy(1) = tmp_iom_output(1)
2198             noos_iom_dummy(2) = tmp_iom_output(2)
2199             noos_iom_dummy(3) = tmp_iom_output(3)
2200             
2201             !noos_iom_dummy(1) = (zsumclasses( 7)+zsumclasses( 8))
2202             !noos_iom_dummy(2) = zsumclasses( 7)
2203             !noos_iom_dummy(3) = zsumclasses( 8)
2204             
2205             if ( lwp .AND. verbose ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
2206             CALL iom_put(noos_var_sect_name,  noos_iom_dummy(:) ) 
2207         ENDIF
2208
2209         noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_salt'
2210         IF (iom_use(noos_var_sect_name)) THEN
2211             noos_iom_dummy(:) = 0.
2212             tmp_iom_output(:) = 0.
2213             
2214             tmp_iom_output(1) = (zsumclasses( 9)+zsumclasses( 10))
2215             tmp_iom_output(2) = zsumclasses( 9)
2216             tmp_iom_output(3) = zsumclasses( 10)
2217             
2218             ! Convert to  MT/s
2219             tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
2220             tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
2221             tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
2222             
2223             
2224             ! limit maximum and minimum values in iom_put
2225             if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
2226             if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
2227             if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
2228             if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
2229             if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
2230             if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
2231             
2232             ! Set NaN's to Zero         
2233             if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
2234             if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
2235             if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
2236             
2237             noos_iom_dummy(1) = tmp_iom_output(1)
2238             noos_iom_dummy(2) = tmp_iom_output(2)
2239             noos_iom_dummy(3) = tmp_iom_output(3)
2240             
2241             !noos_iom_dummy(1) = (zsumclasses( 9)+zsumclasses( 10))
2242             !noos_iom_dummy(2) = zsumclasses( 9)
2243             !noos_iom_dummy(3) = zsumclasses( 10)
2244             
2245             if ( lwp .AND. verbose ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
2246             CALL iom_put(noos_var_sect_name,  noos_iom_dummy(:) )
2247             noos_iom_dummy(:) = 0.         
2248             tmp_iom_output(:) = 0.
2249        ENDIF
2250
2251        DEALLOCATE(noos_iom_dummy)
2252     ENDIF
2253     
2254
2255     DO jclass=1,MAX(1,sec%nb_class-1)
2256   
2257        classe   = 'N       '
2258        zbnd1   = 0._wp
2259        zbnd2   = 0._wp
2260
2261        !insitu density classes transports
2262        IF( ( sec%zsigi(jclass)   .NE. 99._wp ) .AND. &
2263            ( sec%zsigi(jclass+1) .NE. 99._wp )       )THEN
2264           classe = 'DI       '
2265           zbnd1 = sec%zsigi(jclass)
2266           zbnd2 = sec%zsigi(jclass+1)
2267        ENDIF
2268        !potential density classes transports
2269        IF( ( sec%zsigp(jclass)   .NE. 99._wp ) .AND. &
2270            ( sec%zsigp(jclass+1) .NE. 99._wp )       )THEN
2271           classe = 'DP      '
2272           zbnd1 = sec%zsigp(jclass)
2273           zbnd2 = sec%zsigp(jclass+1)
2274        ENDIF
2275        !depth classes transports
2276        IF( ( sec%zlay(jclass)    .NE. 99._wp ) .AND. &
2277            ( sec%zlay(jclass+1)  .NE. 99._wp )       )THEN
2278           classe = 'Z       '
2279           zbnd1 = sec%zlay(jclass)
2280           zbnd2 = sec%zlay(jclass+1)
2281        ENDIF
2282        !salinity classes transports
2283        IF( ( sec%zsal(jclass) .NE. 99._wp    ) .AND. &
2284            ( sec%zsal(jclass+1) .NE. 99._wp  )       )THEN
2285           classe = 'S       '
2286           zbnd1 = sec%zsal(jclass)
2287           zbnd2 = sec%zsal(jclass+1)   
2288        ENDIF
2289        !temperature classes transports
2290        IF( ( sec%ztem(jclass) .NE. 99._wp     ) .AND. &
2291            ( sec%ztem(jclass+1) .NE. 99._wp     )       ) THEN
2292           classe = 'T       '
2293           zbnd1 = sec%ztem(jclass)
2294           zbnd2 = sec%ztem(jclass+1)
2295        ENDIF
2296                 
2297        !write volume transport per class
2298        IF( lwp ) THEN
2299           
2300            IF  ( ln_dct_ascii ) THEN
2301                CALL FLUSH(numdct_NOOS)
2302
2303                !WRITE(numdct_NOOS,'(9e12.4E2)')   sec%transport( 1,jclass)+sec%transport( 2,jclass) , sec%transport( 1,jclass), sec%transport( 2,jclass), &
2304                !                                  sec%transport( 7,jclass)+sec%transport( 8,jclass) , sec%transport( 7,jclass), sec%transport( 8,jclass), &
2305                !                                  sec%transport( 9,jclass)+sec%transport(10,jclass) , sec%transport( 9,jclass), sec%transport(10,jclass)
2306                WRITE(numdct_NOOS,'(3F18.3,6e16.8E2)')   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            ELSE
2310
2311                CALL FLUSH(numdct_NOOS)
2312                WRITE(numdct_NOOS)   sec%transport( 1,jclass)+sec%transport( 2,jclass) , sec%transport( 1,jclass), sec%transport( 2,jclass), &
2313                                     sec%transport( 7,jclass)+sec%transport( 8,jclass) , sec%transport( 7,jclass), sec%transport( 8,jclass), &
2314                                     sec%transport( 9,jclass)+sec%transport(10,jclass) , sec%transport( 9,jclass), sec%transport(10,jclass)
2315            ENDIF
2316        ENDIF
2317
2318     ENDDO
2319     
2320     !IF  ( ln_dct_ascii ) THEN
2321        if ( lwp ) CALL FLUSH(numdct_NOOS)
2322     !ENDIF
2323
2324     !CALL wrk_dealloc(nb_type_class , zsumclasses ) 
2325     DEALLOCATE(  zsumclasses )
2326
2327  END SUBROUTINE dia_dct_wri_NOOS
2328
2329
2330
2331
2332
2333
2334
2335  SUBROUTINE dia_dct_wri_NOOS_iom(kt,ksec,sec)
2336     !!-------------------------------------------------------------
2337     !! Write transport output in numdct using NOOS formatting
2338     !!
2339     !! Purpose: Write  transports in ascii files
2340     !!
2341     !! Method:
2342     !!        1. Write volume transports in "volume_transport"
2343     !!           Unit: Sv : area * Velocity / 1.e6
2344     !!
2345     !!        2. Write heat transports in "heat_transport"
2346     !!           Unit: Peta W : area * Velocity * T * rhau * Cp / 1.e15
2347     !!
2348     !!        3. Write salt transports in "salt_transport"
2349     !!           Unit: 10^9 g m^3 / s : area * Velocity * S / 1.e6
2350     !!
2351     !!-------------------------------------------------------------
2352     !!arguments
2353     INTEGER, INTENT(IN)          :: kt          ! time-step
2354     TYPE(SECTION), INTENT(INOUT) :: sec         ! section to write   
2355     INTEGER ,INTENT(IN)          :: ksec        ! section number
2356
2357     !!local declarations
2358     INTEGER               :: jclass,ji             ! Dummy loop
2359     CHARACTER(len=2)      :: classe             ! Classname
2360     REAL(wp)              :: zbnd1,zbnd2        ! Class bounds
2361     !REAL(wp)              :: zslope             ! section's slope coeff
2362     !
2363     REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace
2364     CHARACTER(len=3)      :: noos_sect_name             ! Classname
2365     CHARACTER(len=25)      :: noos_var_sect_name
2366     REAL(wp), ALLOCATABLE, DIMENSION(:) ::   noos_iom_dummy
2367     INTEGER               :: IERR
2368     
2369     REAL(wp), DIMENSION(3) :: tmp_iom_output
2370     REAL(wp)               :: max_iom_val
2371     LOGICAL       ::   verbose     
2372     verbose = ln_dct_verbose! .false.
2373     
2374     !!-------------------------------------------------------------
2375     
2376     
2377     
2378     IF( lwp .AND. verbose ) THEN
2379        WRITE(numout,*) " "
2380        WRITE(numout,*) "dia_dct_wri_NOOS_IOM: write transports through sections at timestep: ", kt
2381        WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
2382     ENDIF   
2383       
2384     !CALL wrk_alloc(nb_type_class , zsumclasses ) 
2385     ALLOCATE( zsumclasses(nb_type_class),  STAT= ierr )
2386     IF( ierr /= 0 )   CALL ctl_stop( 'dia_dct_wri_NOOS_IOM: failed to allocate zsumclasses array' )
2387
2388     zsumclasses(:)=0._wp
2389
2390!   
2391!     ! Sum all classes together, to give one values per type (pos tran, neg vol trans etc...).
2392!     DO jclass=1,MAX(1,sec%nb_class-1)
2393!        zsumclasses(1:nb_type_class)=zsumclasses(1:nb_type_class)+sec%transport(1:nb_type_class,jclass)
2394!     ENDDO
2395!
2396!     classe   = 'total   '
2397     zbnd1   = 0._wp
2398     zbnd2   = 0._wp
2399     
2400     zsumclasses(1) = transports_3d_inst_sum(1,ksec,2)
2401     zsumclasses(2) = transports_3d_inst_sum(1,ksec,3)
2402     zsumclasses(7) = transports_3d_inst_sum(4,ksec,2)
2403     zsumclasses(8) = transports_3d_inst_sum(4,ksec,3)
2404     zsumclasses(9) = transports_3d_inst_sum(5,ksec,2)
2405     zsumclasses(10) = transports_3d_inst_sum(5,ksec,3)
2406
2407     
2408     write (noos_sect_name, "(I0.3)")  ksec
2409     
2410     !IF ( nn_dct_iom_cont .EQ. 2 ) THEN
2411     max_iom_val = 1.e10
2412     ALLOCATE( noos_iom_dummy(3),  STAT= ierr )
2413        IF( ierr /= 0 )   CALL ctl_stop( 'dia_dct_wri_NOOS: failed to allocate noos_iom_dummy array' )
2414
2415     noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_trans'
2416     IF (iom_use(noos_var_sect_name)) THEN
2417         noos_iom_dummy(:) = 0.
2418         tmp_iom_output(:) = 0.
2419         
2420         tmp_iom_output(1) = (zsumclasses( 1)+zsumclasses( 2))
2421         tmp_iom_output(2) = zsumclasses( 1)
2422         tmp_iom_output(3) = zsumclasses( 2)
2423         
2424         ! Convert to Sv
2425         tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
2426         tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
2427         tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
2428         
2429         ! limit maximum and minimum values in iom_put
2430         if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
2431         if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
2432         if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
2433         if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
2434         if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
2435         if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
2436         
2437         ! Set NaN's to Zero         
2438         if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
2439         if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
2440         if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
2441         
2442         noos_iom_dummy(1) = tmp_iom_output(1)
2443         noos_iom_dummy(2) = tmp_iom_output(2)
2444         noos_iom_dummy(3) = tmp_iom_output(3)
2445         
2446         !noos_iom_dummy(1) = (zsumclasses( 1)+zsumclasses( 2))
2447         !noos_iom_dummy(2) = zsumclasses( 1)
2448         !noos_iom_dummy(3) = zsumclasses( 2)
2449         
2450         
2451         
2452         if ( lwp .AND. verbose ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name,tmp_iom_output(1)
2453         CALL iom_put( noos_var_sect_name,  noos_iom_dummy(:)  )
2454     ENDIF
2455
2456     noos_var_sect_name =  "noos_" // trim(noos_sect_name) // '_heat'
2457     IF (iom_use(noos_var_sect_name)) THEN
2458         noos_iom_dummy(:) = 0.
2459         tmp_iom_output(:) = 0.
2460         
2461         tmp_iom_output(1) = (zsumclasses( 7)+zsumclasses( 8))
2462         tmp_iom_output(2) = zsumclasses( 7)
2463         tmp_iom_output(3) = zsumclasses( 8)
2464         
2465         ! Convert to TJ/s
2466         tmp_iom_output(1) = tmp_iom_output(1)*1.E-12
2467         tmp_iom_output(2) = tmp_iom_output(2)*1.E-12
2468         tmp_iom_output(3) = tmp_iom_output(3)*1.E-12
2469         
2470         ! limit maximum and minimum values in iom_put
2471         if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
2472         if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
2473         if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
2474         if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
2475         if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
2476         if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
2477         
2478         ! Set NaN's to Zero         
2479         if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
2480         if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
2481         if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
2482         
2483         noos_iom_dummy(1) = tmp_iom_output(1)
2484         noos_iom_dummy(2) = tmp_iom_output(2)
2485         noos_iom_dummy(3) = tmp_iom_output(3)
2486         
2487         !noos_iom_dummy(1) = (zsumclasses( 7)+zsumclasses( 8))
2488         !noos_iom_dummy(2) = zsumclasses( 7)
2489         !noos_iom_dummy(3) = zsumclasses( 8)
2490         
2491         if ( lwp .AND. verbose ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name,tmp_iom_output(1)
2492         CALL iom_put(noos_var_sect_name,  noos_iom_dummy(:) ) 
2493     ENDIF
2494
2495     noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_salt'
2496     IF (iom_use(noos_var_sect_name)) THEN
2497         noos_iom_dummy(:) = 0.
2498         tmp_iom_output(:) = 0.
2499         
2500         tmp_iom_output(1) = (zsumclasses( 9)+zsumclasses( 10))
2501         tmp_iom_output(2) = zsumclasses( 9)
2502         tmp_iom_output(3) = zsumclasses( 10)
2503         
2504         ! Convert to  MT/s
2505         tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
2506         tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
2507         tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
2508         
2509         
2510         ! limit maximum and minimum values in iom_put
2511         if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
2512         if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
2513         if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
2514         if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
2515         if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
2516         if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
2517         
2518         ! Set NaN's to Zero         
2519         if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
2520         if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
2521         if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
2522         
2523         noos_iom_dummy(1) = tmp_iom_output(1)
2524         noos_iom_dummy(2) = tmp_iom_output(2)
2525         noos_iom_dummy(3) = tmp_iom_output(3)
2526         
2527         !noos_iom_dummy(1) = (zsumclasses( 9)+zsumclasses( 10))
2528         !noos_iom_dummy(2) = zsumclasses( 9)
2529         !noos_iom_dummy(3) = zsumclasses( 10)
2530         
2531         if ( lwp .AND. verbose ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name,tmp_iom_output(1)
2532         CALL iom_put(noos_var_sect_name,  noos_iom_dummy(:)  )
2533         noos_iom_dummy(:) = 0.         
2534         tmp_iom_output(:) = 0.
2535    ENDIF
2536
2537    DEALLOCATE(noos_iom_dummy)
2538     !ENDIF
2539     
2540     
2541     !CALL wrk_dealloc(nb_type_class , zsumclasses ) 
2542     DEALLOCATE(  zsumclasses )
2543
2544  END SUBROUTINE dia_dct_wri_NOOS_iom
2545
2546
2547  SUBROUTINE dia_dct_wri_NOOS_h(hr,ksec,sec)
2548     !!-------------------------------------------------------------
2549     !! As routine dia_dct_wri_NOOS but for hourly output files
2550     !!
2551     !! Write transport output in numdct using NOOS formatting
2552     !!
2553     !! Purpose: Write  transports in ascii files
2554     !!
2555     !! Method:
2556     !!        1. Write volume transports in "volume_transport"
2557     !!           Unit: Sv : area * Velocity / 1.e6
2558     !!
2559     !!-------------------------------------------------------------
2560     !!arguments
2561     INTEGER, INTENT(IN)          :: hr          ! hour => effectively kt/12
2562     TYPE(SECTION), INTENT(INOUT) :: sec         ! section to write   
2563     INTEGER ,INTENT(IN)          :: ksec        ! section number
2564
2565     !!local declarations
2566     INTEGER               :: jclass,jhr            ! Dummy loop
2567     CHARACTER(len=2)      :: classe             ! Classname
2568     REAL(wp)              :: zbnd1,zbnd2        ! Class bounds
2569     REAL(wp)              :: zslope             ! section's slope coeff
2570     !
2571     REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace
2572     CHARACTER(len=3)      :: noos_sect_name             ! Classname
2573     CHARACTER(len=25)      :: noos_var_sect_name
2574     REAL(wp), ALLOCATABLE, DIMENSION(:) ::   noos_iom_dummy
2575     INTEGER               :: IERR
2576     LOGICAL       ::   verbose     
2577     verbose = ln_dct_verbose! .false.
2578     
2579     !!-------------------------------------------------------------
2580
2581     IF( lwp .AND. verbose ) THEN
2582        WRITE(numout,*) " "
2583        WRITE(numout,*) "dia_dct_wri_NOOS_h: write transports through section Transect:",ksec-1," at timestep: ", hr
2584        WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
2585     ENDIF   
2586     
2587     !CALL wrk_alloc(nb_type_class , zsumclasses )
2588     ALLOCATE( zsumclasses(nb_type_class),  STAT= ierr )
2589     IF( ierr /= 0 )   CALL ctl_stop( 'dia_dct_wri_NOOS_h: failed to allocate zsumclasses array' )
2590     
2591     
2592     write (noos_sect_name, "(I03)")  ksec
2593     
2594     ALLOCATE( noos_iom_dummy(3),  STAT= ierr )
2595        IF( ierr /= 0 )   CALL ctl_stop( 'dia_dct_wri_NOOS_h: failed to allocate noos_iom_dummy array' )
2596
2597
2598
2599     
2600
2601     zsumclasses(:)=0._wp
2602     zslope = sec%slopeSection       
2603
2604     ! Sum up all classes, to give the total per type (pos vol trans, neg vol trans etc...)
2605     DO jclass=1,MAX(1,sec%nb_class-1)
2606        zsumclasses(1:nb_type_class)=zsumclasses(1:nb_type_class)+sec%transport_h(1:nb_type_class,jclass)
2607     ENDDO
2608       
2609     
2610     ! JT I think changing the sign of output according to the zslope is redundant
2611     
2612     !write volume transport per class
2613     ! Sum positive and vol trans for all classes in first cell of array
2614
2615     z_hr_output(ksec,1,1)= (zsumclasses(1)+zsumclasses(2))
2616     z_hr_output(ksec,2,1)= zsumclasses(1)
2617     z_hr_output(ksec,3,1)= zsumclasses(2)
2618
2619     ! Sum positive and vol trans for each classes in following cell of array
2620     DO jclass=1,MAX(1,sec%nb_class-1)
2621        z_hr_output(ksec,1,jclass+1)= (sec%transport_h(1,jclass)+sec%transport_h(2,jclass))
2622        z_hr_output(ksec,2,jclass+1)= sec%transport_h(1,jclass)
2623        z_hr_output(ksec,3,jclass+1)= sec%transport_h(2,jclass)
2624     ENDDO
2625
2626   
2627    IF( lwp )  THEN
2628        ! JT IF ( hr .eq. 48._wp ) THEN
2629        ! 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
2630        ! JT    DO jhr=25,48
2631        ! 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) )
2632        ! JT    ENDDO
2633        ! JT ENDIF
2634
2635
2636
2637        IF ( ln_dct_ascii ) THEN
2638            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
2639            WRITE(numdct_NOOS_h,'(11F18.3)')  z_hr_output(ksec,1,1), (z_hr_output(ksec,1,jclass+1), jclass=1,MAX(1,10) )
2640            WRITE(numdct_NOOS_h,'(11F18.3)')  z_hr_output(ksec,2,1), (z_hr_output(ksec,2,jclass+1), jclass=1,MAX(1,10) )
2641            WRITE(numdct_NOOS_h,'(11F18.3)')  z_hr_output(ksec,3,1), (z_hr_output(ksec,3,jclass+1), jclass=1,MAX(1,10) )
2642            CALL FLUSH(numdct_NOOS_h)
2643        ELSE
2644            WRITE(numdct_NOOS_h) nyear,nmonth,nday,MOD(hr,24),ksec-1,sec%nb_class-1
2645            WRITE(numdct_NOOS_h)  z_hr_output(ksec,1,1), (z_hr_output(ksec,1,jclass+1), jclass=1,MAX(1,10) )
2646            WRITE(numdct_NOOS_h)  z_hr_output(ksec,2,1), (z_hr_output(ksec,2,jclass+1), jclass=1,MAX(1,10) )
2647            WRITE(numdct_NOOS_h)  z_hr_output(ksec,3,1), (z_hr_output(ksec,3,jclass+1), jclass=1,MAX(1,10) )
2648            CALL FLUSH(numdct_NOOS_h)
2649        ENDIF
2650
2651
2652     ENDIF 
2653
2654
2655     !CALL wrk_dealloc(nb_type_class , zsumclasses )
2656     DEALLOCATE( zsumclasses )
2657     
2658     DEALLOCATE(noos_iom_dummy)
2659
2660
2661
2662  END SUBROUTINE dia_dct_wri_NOOS_h
2663
2664  SUBROUTINE dia_dct_wri(kt,ksec,sec)
2665     !!-------------------------------------------------------------
2666     !! Write transport output in numdct
2667     !!
2668     !! Purpose: Write  transports in ascii files
2669     !!
2670     !! Method:
2671     !!        1. Write volume transports in "volume_transport"
2672     !!           Unit: Sv : area * Velocity / 1.e6
2673     !!
2674     !!        2. Write heat transports in "heat_transport"
2675     !!           Unit: Peta W : area * Velocity * T * rhau * Cp / 1.e15
2676     !!
2677     !!        3. Write salt transports in "salt_transport"
2678     !!           Unit: 10^9 g m^3 / s : area * Velocity * S / 1.e6
2679     !!
2680     !!-------------------------------------------------------------
2681     !!arguments
2682     INTEGER, INTENT(IN)          :: kt          ! time-step
2683     TYPE(SECTION), INTENT(INOUT) :: sec         ! section to write   
2684     INTEGER ,INTENT(IN)          :: ksec        ! section number
2685
2686     !!local declarations
2687     INTEGER                            :: ierr  ! error for allocate
2688     INTEGER               :: jclass             ! Dummy loop
2689     CHARACTER(len=2)      :: classe             ! Classname
2690     REAL(wp)              :: zbnd1,zbnd2        ! Class bounds
2691     REAL(wp)              :: zslope             ! section's slope coeff
2692     !
2693     REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace
2694     !!-------------------------------------------------------------
2695     !CALL wrk_alloc(nb_type_class , zsumclasses ) 
2696     ALLOCATE( zsumclasses(nb_type_class),  STAT= ierr )
2697     IF( ierr /= 0 )   CALL ctl_stop( 'dia_dct_wri: failed to allocate zsumclasses array' )
2698
2699     zsumclasses(:)=0._wp
2700     zslope = sec%slopeSection       
2701
2702 
2703     DO jclass=1,MAX(1,sec%nb_class-1)
2704
2705        classe   = 'N       '
2706        zbnd1   = 0._wp
2707        zbnd2   = 0._wp
2708        zsumclasses(1:nb_type_class)=zsumclasses(1:nb_type_class)+sec%transport(1:nb_type_class,jclass)
2709
2710   
2711        !insitu density classes transports
2712        IF( ( sec%zsigi(jclass)   .NE. 99._wp ) .AND. &
2713            ( sec%zsigi(jclass+1) .NE. 99._wp )       )THEN
2714           classe = 'DI       '
2715           zbnd1 = sec%zsigi(jclass)
2716           zbnd2 = sec%zsigi(jclass+1)
2717        ENDIF
2718        !potential density classes transports
2719        IF( ( sec%zsigp(jclass)   .NE. 99._wp ) .AND. &
2720            ( sec%zsigp(jclass+1) .NE. 99._wp )       )THEN
2721           classe = 'DP      '
2722           zbnd1 = sec%zsigp(jclass)
2723           zbnd2 = sec%zsigp(jclass+1)
2724        ENDIF
2725        !depth classes transports
2726        IF( ( sec%zlay(jclass)    .NE. 99._wp ) .AND. &
2727            ( sec%zlay(jclass+1)  .NE. 99._wp )       )THEN
2728           classe = 'Z       '
2729           zbnd1 = sec%zlay(jclass)
2730           zbnd2 = sec%zlay(jclass+1)
2731        ENDIF
2732        !salinity classes transports
2733        IF( ( sec%zsal(jclass) .NE. 99._wp    ) .AND. &
2734            ( sec%zsal(jclass+1) .NE. 99._wp  )       )THEN
2735           classe = 'S       '
2736           zbnd1 = sec%zsal(jclass)
2737           zbnd2 = sec%zsal(jclass+1)   
2738        ENDIF
2739        !temperature classes transports
2740        IF( ( sec%ztem(jclass) .NE. 99._wp     ) .AND. &
2741            ( sec%ztem(jclass+1) .NE. 99._wp     )       ) THEN
2742           classe = 'T       '
2743           zbnd1 = sec%ztem(jclass)
2744           zbnd2 = sec%ztem(jclass+1)
2745        ENDIF
2746                 
2747        !write volume transport per class
2748        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, &
2749                              jclass,classe,zbnd1,zbnd2,&
2750                              sec%transport(1,jclass),sec%transport(2,jclass), &
2751                              sec%transport(1,jclass)+sec%transport(2,jclass)
2752
2753        IF( sec%llstrpond )THEN
2754
2755           !write heat transport per class:
2756           WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope,  &
2757                              jclass,classe,zbnd1,zbnd2,&
2758                              sec%transport(7,jclass)*1000._wp*rcp/1.e15,sec%transport(8,jclass)*1000._wp*rcp/1.e15, &
2759                              ( sec%transport(7,jclass)+sec%transport(8,jclass) )*1000._wp*rcp/1.e15
2760           !write salt transport per class
2761           WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope,  &
2762                              jclass,classe,zbnd1,zbnd2,&
2763                              sec%transport(9,jclass)*1000._wp/1.e9,sec%transport(10,jclass)*1000._wp/1.e9,&
2764                              (sec%transport(9,jclass)+sec%transport(10,jclass))*1000._wp/1.e9
2765        ENDIF
2766
2767     ENDDO
2768
2769     zbnd1 = 0._wp
2770     zbnd2 = 0._wp
2771     jclass=0
2772
2773     !write total volume transport
2774     WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, &
2775                           jclass,"total",zbnd1,zbnd2,&
2776                           zsumclasses(1),zsumclasses(2),zsumclasses(1)+zsumclasses(2)
2777
2778     IF( sec%llstrpond )THEN
2779
2780        !write total heat transport
2781        WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope, &
2782                           jclass,"total",zbnd1,zbnd2,&
2783                           zsumclasses(7)* 1000._wp*rcp/1.e15,zsumclasses(8)* 1000._wp*rcp/1.e15,&
2784                           (zsumclasses(7)+zsumclasses(8) )* 1000._wp*rcp/1.e15
2785        !write total salt transport
2786        WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope, &
2787                           jclass,"total",zbnd1,zbnd2,&
2788                           zsumclasses(9)*1000._wp/1.e9,zsumclasses(10)*1000._wp/1.e9,&
2789                           (zsumclasses(9)+zsumclasses(10))*1000._wp/1.e9
2790     ENDIF
2791
2792     
2793     IF ( sec%ll_ice_section) THEN
2794        !write total ice volume transport
2795        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,&
2796                              jclass,"ice_vol",zbnd1,zbnd2,&
2797                              sec%transport(11,1),sec%transport(12,1),&
2798                              sec%transport(11,1)+sec%transport(12,1)
2799        !write total ice surface transport
2800        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,&
2801                              jclass,"ice_surf",zbnd1,zbnd2,&
2802                              sec%transport(13,1),sec%transport(14,1), &
2803                              sec%transport(13,1)+sec%transport(14,1) 
2804     ENDIF
2805                                             
2806118 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3F12.4)
2807119 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6)
2808
2809     !CALL wrk_dealloc(nb_type_class , zsumclasses ) 
2810     DEALLOCATE ( zsumclasses ) 
2811  END SUBROUTINE dia_dct_wri
2812
2813   PURE  FUNCTION interp(ki, kj, kk, cd_point,var) 
2814  !!----------------------------------------------------------------------
2815  !!
2816  !!   Purpose: compute temperature/salinity/density at U-point or V-point
2817  !!   --------
2818  !!
2819  !!   Method:
2820  !!   ------
2821  !!
2822  !!   ====> full step and partial step
2823  !!
2824  !!
2825  !!    |    I          |    I+1           |    Z=temperature/salinity/density at U-poinT
2826  !!    |               |                  |
2827  !!  ----------------------------------------  1. Veritcale interpolation: compute zbis
2828  !!    |               |                  |       interpolation between ptab(I,J,K) and ptab(I,J,K+1)
2829  !!    |               |                  |       zbis =
2830  !!    |               |                  |      [ e3w(I+1,J,K)*ptab(I,J,K) + ( e3w(I,J,K) - e3w(I+1,J,K) ) * ptab(I,J,K-1) ]
2831  !!    |               |                  |      /[ e3w(I+1,J,K) + e3w(I,J,K) - e3w(I+1,J,K) ]
2832  !!    |               |                  |
2833  !!    |               |                  |    2. Horizontal interpolation: compute value at U/V point
2834  !!K-1 | ptab(I,J,K-1) |                  |       interpolation between zbis and ptab(I+1,J,K) 
2835  !!    |     .         |                  |
2836  !!    |     .         |                  |       interp = ( 0.5*zet2*zbis + 0.5*zet1*ptab(I+1,J,K) )/(0.5*zet2+0.5*zet1)
2837  !!    |     .         |                  |
2838  !!  ------------------------------------------
2839  !!    |     .         |                  |
2840  !!    |     .         |                  |
2841  !!    |     .         |                  |
2842  !!K   |    zbis.......U...ptab(I+1,J,K)  |
2843  !!    |     .         |                  |
2844  !!    | ptab(I,J,K)   |                  |
2845  !!    |               |------------------|
2846  !!    |               | partials         |
2847  !!    |               |  steps           |
2848  !!  -------------------------------------------
2849  !!    <----zet1------><----zet2--------->
2850  !!
2851  !!
2852  !!   ====>  s-coordinate
2853  !!     
2854  !!    |                |                  |   1. Compute distance between T1 and U points: SQRT( zdep1^2 + (0.5 * zet1 )^2
2855  !!    |                |                  |      Compute distance between T2 and U points: SQRT( zdep2^2 + (0.5 * zet2 )^2
2856  !!    |                | ptab(I+1,J,K)    |
2857  !!    |                |      T2          |   2. Interpolation between  T1 and T2 values at U point
2858  !!    |                |      ^           |   
2859  !!    |                |      | zdep2     |   
2860  !!    |                |      |           |   
2861  !!    |       ^        U      v           |
2862  !!    |       |        |                  |
2863  !!    |       | zdep1  |                  |   
2864  !!    |       v        |                  |
2865  !!    |      T1        |                  |
2866  !!    | ptab(I,J,K)    |                  |
2867  !!    |                |                  |
2868  !!    |                |                  |
2869  !!
2870  !!    <----zet1--------><----zet2--------->
2871  !!
2872  !!----------------------------------------------------------------------
2873  !*arguments
2874  INTEGER, INTENT(IN)                          :: ki, kj, kk   ! coordinate of point
2875  INTEGER, INTENT(IN)                          :: var   !  which variable
2876  CHARACTER(len=1), INTENT(IN)                 :: cd_point     ! type of point (U, V)
2877  REAL(wp)                                     :: interp       ! interpolated variable
2878
2879  !*local declations
2880  INTEGER :: ii1, ij1, ii2, ij2                                ! local integer
2881  REAL(wp):: ze3t, zfse3, zwgt1, zwgt2, zbis, zdepu            ! local real
2882  REAL(wp):: zet1, zet2                                        ! weight for interpolation
2883  REAL(wp):: zdep1,zdep2                                       ! differences of depth
2884  REAL(wp):: zmsk                                              ! mask value
2885  !!----------------------------------------------------------------------
2886
2887 
2888
2889  IF( cd_point=='U' )THEN
2890     ii1 = ki    ; ij1 = kj 
2891     ii2 = ki+1  ; ij2 = kj 
2892
2893     zet1=e1t(ii1,ij1)
2894     zet2=e1t(ii2,ij2)
2895     zmsk=umask(ii1,ij1,kk)
2896 
2897
2898  ELSE ! cd_point=='V'
2899     ii1 = ki    ; ij1 = kj 
2900     ii2 = ki    ; ij2 = kj+1 
2901
2902     zet1=e2t(ii1,ij1)
2903     zet2=e2t(ii2,ij2)
2904     zmsk=vmask(ii1,ij1,kk)
2905
2906  ENDIF
2907
2908  IF( ln_sco )THEN   ! s-coordinate case
2909
2910     zdepu = ( gdept_n(ii1,ij1,kk) +  gdept_n(ii2,ij2,kk) ) /2 
2911     zdep1 = gdept_n(ii1,ij1,kk) - zdepu
2912     zdep2 = gdept_n(ii2,ij2,kk) - zdepu
2913
2914     ! weights
2915     zwgt1 = SQRT( ( 0.5 * zet1 ) * ( 0.5 * zet1 ) + ( zdep1 * zdep1 ) )
2916     zwgt2 = SQRT( ( 0.5 * zet2 ) * ( 0.5 * zet2 ) + ( zdep2 * zdep2 ) )
2917 
2918     ! result
2919           SELECT CASE( var )
2920              CASE(0)  ;    interp = zmsk * ( zwgt2 *  tsn(ii1,ij1,kk,jp_tem) + zwgt1 *  tsn(ii1,ij1,kk,jp_tem) ) / ( zwgt2 + zwgt1 )
2921              CASE(1)  ;    interp = zmsk * ( zwgt2 *  tsn(ii1,ij1,kk,jp_sal) + zwgt1 *  tsn(ii1,ij1,kk,jp_sal) ) / ( zwgt2 + zwgt1 )
2922              CASE(2)  ;    interp = zmsk * ( zwgt2 *  rhop(ii1,ij1,kk) + zwgt1 *  rhop(ii1,ij1,kk) ) / ( zwgt2 + zwgt1 )
2923              CASE(3)  ;    interp = zmsk * ( zwgt2 *  (rhd(ii1,ij1,kk)*rau0+rau0) + zwgt1 *  (rhd(ii1,ij1,kk)*rau0+rau0) ) / ( zwgt2 + zwgt1 )
2924           END SELECT
2925
2926  ELSE       ! full step or partial step case
2927
2928#if defined key_vvl
2929
2930     !ze3t  = fse3t_n(ii2,ij2,kk) - fse3t_n(ii1,ij1,kk)
2931     !zwgt1 = ( fse3w_n(ii2,ij2,kk) - fse3w_n(ii1,ij1,kk) ) / fse3w_n(ii2,ij2,kk)
2932     !zwgt2 = ( fse3w_n(ii1,ij1,kk) - fse3w_n(ii2,ij2,kk) ) / fse3w_n(ii1,ij1,kk)
2933
2934     ze3t  = e3t_n(ii2,ij2,kk) - e3t_n(ii1,ij1,kk) 
2935     zwgt1 = ( e3w_n(ii2,ij2,kk) - e3w_n(ii1,ij1,kk) ) / e3w_n(ii2,ij2,kk)
2936     zwgt2 = ( e3w_n(ii1,ij1,kk) - e3w_n(ii2,ij2,kk) ) / e3w_n(ii1,ij1,kk)
2937
2938#else
2939
2940     !ze3t  = fse3t(ii2,ij2,kk)   - fse3t(ii1,ij1,kk)
2941     !zwgt1 = ( fse3w(ii2,ij2,kk) - fse3w(ii1,ij1,kk) ) / fse3w(ii2,ij2,kk)
2942     !zwgt2 = ( fse3w(ii1,ij1,kk) - fse3w(ii2,ij2,kk) ) / fse3w(ii1,ij1,kk)
2943
2944     !ze3t  = e3t(ii2,ij2,kk)   - e3t(ii1,ij1,kk)
2945     !zwgt1 = ( e3w(ii2,ij2,kk) - e3w(ii1,ij1,kk) ) / e3w(ii2,ij2,kk)
2946     !zwgt2 = ( e3w(ii1,ij1,kk) - e3w(ii2,ij2,kk) ) / e3w(ii1,ij1,kk)
2947
2948
2949     ze3t  = e3t_n(ii2,ij2,kk) - e3t_n(ii1,ij1,kk) 
2950     zwgt1 = ( e3w_n(ii2,ij2,kk) - e3w_n(ii1,ij1,kk) ) / e3w_n(ii2,ij2,kk)
2951     zwgt2 = ( e3w_n(ii1,ij1,kk) - e3w_n(ii2,ij2,kk) ) / e3w_n(ii1,ij1,kk)
2952
2953#endif
2954
2955     IF(kk .NE. 1)THEN
2956
2957        IF( ze3t >= 0. )THEN 
2958           ! zbis
2959           SELECT CASE( var )
2960           CASE(0) 
2961                     zbis   =  tsn(ii2,ij2,kk,jp_tem) + zwgt1 *  (tsn(ii2,ij2,kk-1,jp_tem)-tsn(ii2,ij2,kk,jp_tem)   )
2962                     interp =  zmsk * ( zet2 * tsn(ii1,ij1,kk,jp_tem) + zet1 * zbis )/( zet1 + zet2 )
2963           CASE(1) 
2964                     zbis   =  tsn(ii2,ij2,kk,jp_sal) + zwgt1 *  (tsn(ii2,ij2,kk-1,jp_sal)-tsn(ii2,ij2,kk,jp_sal)   )
2965                     interp =  zmsk * ( zet2 * tsn(ii1,ij1,kk,jp_sal) + zet1 * zbis )/( zet1 + zet2 )
2966           CASE(2) 
2967                     zbis   =  rhop(ii2,ij2,kk) + zwgt1 *  (rhop(ii2,ij2,kk-1)-rhop(ii2,ij2,kk)   )
2968                     interp =  zmsk * ( zet2 * rhop(ii1,ij1,kk) + zet1 * zbis )/( zet1 + zet2 )
2969           CASE(3) 
2970                     zbis   =  (rhd(ii2,ij2,kk)*rau0+rau0) + zwgt1 *  ( (rhd(ii2,ij2,kk-1)*rau0+rau0)-(rhd(ii2,ij2,kk)*rau0+rau0)   )
2971                     interp =  zmsk * ( zet2 * (rhd(ii1,ij1,kk)*rau0+rau0) + zet1 * zbis )/( zet1 + zet2 )
2972           END SELECT
2973           ! result
2974        ELSE
2975           ! zbis
2976           SELECT CASE( var )
2977           CASE(0) 
2978                 zbis   = tsn(ii1,ij1,kk,jp_tem) + zwgt2 * ( tsn(ii1,ij1,kk-1,jp_tem) - tsn(ii1,ij2,kk,jp_tem) )
2979                 interp = zmsk * ( zet2 * zbis + zet1 * tsn(ii2,ij2,kk,jp_tem) )/( zet1 + zet2 )
2980           CASE(1) 
2981                 zbis   = tsn(ii1,ij1,kk,jp_sal) + zwgt2 * ( tsn(ii1,ij1,kk-1,jp_sal) - tsn(ii1,ij2,kk,jp_sal) )
2982                 interp = zmsk * ( zet2 * zbis + zet1 * tsn(ii2,ij2,kk,jp_sal) )/( zet1 + zet2 )
2983           CASE(2) 
2984                 zbis   = rhop(ii1,ij1,kk) + zwgt2 * ( rhop(ii1,ij1,kk-1) - rhop(ii1,ij2,kk) )
2985                 interp = zmsk * ( zet2 * zbis + zet1 * rhop(ii2,ij2,kk) )/( zet1 + zet2 )
2986           CASE(3) 
2987                 zbis   = (rhd(ii1,ij1,kk)*rau0+rau0) + zwgt2 * ( (rhd(ii1,ij1,kk-1)*rau0+rau0) - (rhd(ii1,ij2,kk)*rau0+rau0) )
2988                 interp = zmsk * ( zet2 * zbis + zet1 * (rhd(ii2,ij2,kk)*rau0+rau0) )/( zet1 + zet2 )
2989           END SELECT
2990        ENDIF   
2991
2992     ELSE
2993        SELECT CASE( var )
2994        CASE(0) 
2995           interp = zmsk * (  zet2 * tsn(ii1,ij1,kk,jp_tem) + zet1 * tsn(ii2,ij2,kk,jp_tem) )/( zet1 + zet2 )
2996        CASE(1) 
2997           interp = zmsk * (  zet2 * tsn(ii1,ij1,kk,jp_sal) + zet1 * tsn(ii2,ij2,kk,jp_sal) )/( zet1 + zet2 )
2998        CASE(2) 
2999           interp = zmsk * (  zet2 * rhop(ii1,ij1,kk) + zet1 * rhop(ii2,ij2,kk) )/( zet1 + zet2 )
3000        CASE(3) 
3001           interp = zmsk * (  zet2 * (rhd(ii1,ij1,kk)*rau0+rau0) + zet1 * (rhd(ii2,ij2,kk)*rau0+rau0) )/( zet1 + zet2 )
3002        END SELECT
3003     ENDIF
3004
3005  ENDIF
3006
3007  END FUNCTION interp
3008
3009!#else
3010!   !!----------------------------------------------------------------------
3011!   !!   Default option :                                       Dummy module
3012!   !!----------------------------------------------------------------------
3013!   LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .FALSE.    !: diamht flag
3014!   PUBLIC
3015!   !! $Id$
3016!CONTAINS
3017
3018!   SUBROUTINE dia_dct_init          ! Dummy routine
3019!      WRITE(*,*) 'dia_dct_init: You should not have seen this print! error?', kt
3020!   END SUBROUTINE dia_dct_init
3021
3022!   SUBROUTINE dia_dct( kt )         ! Dummy routine
3023!      INTEGER, INTENT( in ) :: kt   ! ocean time-step index
3024!      WRITE(*,*) 'dia_dct: You should not have seen this print! error?', kt
3025!   END SUBROUTINE dia_dct
3026!#endif
3027
3028END MODULE diadct
Note: See TracBrowser for help on using the repository browser.