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

Last change on this file since 15378 was 15378, checked in by hadjt, 12 months ago

Updated Region mean and noos cross-sections, but may be slower.
Region mean:
no dummy IOM_put
verbose namelist switch
Region mean and noos consistent heat and salt calc.

NOOS:
instanteous transport arrays added, transport_3d_inst
logical switch in tranport subroutine to update transport_3d, or use transport_3d_inst
New simplified routine to pass instantaneous tranpost to xios via new subroutine

dia_dct_wri_NOOS_iom

Namelist integer to output iom via new routine, old routine, or neither
Verbose name list switch
...

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