source: branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis3/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 11085

Last change on this file since 11085 was 11085, checked in by rrenshaw, 18 months ago

code fixes

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