source: branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 11639

Last change on this file since 11639 was 11639, checked in by rrenshaw, 14 months ago

regional mean and transport diagnostics added for reanalysis

  • Property svn:executable set to *
File size: 132.4 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,*) 'reached end of section_ijglobal.diadct. iost = ',iost, &
619                          ', number of sections read = ', jsec-1
620          EXIT !end of file
621        ENDIF
622       
623       
624        WRITE(cltmp,'(a,i4.4,a,i4.4)')'diadct: read sections : Problem of section number: isec= ',isec,' and jsec= ',jsec
625       
626       
627        IF( jsec .NE. isec )  CALL ctl_stop( cltmp )
628
629        READ(numdct_in)secs(jsec)%name
630        READ(numdct_in)secs(jsec)%llstrpond
631        READ(numdct_in)secs(jsec)%ll_ice_section
632        READ(numdct_in)secs(jsec)%ll_date_line
633        READ(numdct_in)secs(jsec)%coordSec
634        READ(numdct_in)secs(jsec)%nb_class
635        READ(numdct_in)secs(jsec)%zsigi
636        READ(numdct_in)secs(jsec)%zsigp
637        READ(numdct_in)secs(jsec)%zsal
638        READ(numdct_in)secs(jsec)%ztem
639        READ(numdct_in)secs(jsec)%zlay
640        READ(numdct_in)secs(jsec)%slopeSection
641        READ(numdct_in)iptglo
642
643        IF ( ln_NOOS .AND. verbose ) THEN
644           WRITE(numout,*) 'Section name = ',TRIM(secs(jsec)%name)
645           WRITE(numout,*) 'coordSec = ',secs(jsec)%coordSec
646           WRITE(numout,*) 'iptglo = ',iptglo
647        ENDIF
648
649        !debug
650        !-----
651
652        IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
653         
654            WRITE(clformat,'(a,i2,a)') '(A40,', nb_class_max,'(f8.3,1X))' 
655
656            WRITE(numout,*)       "   Section name :                       ",TRIM(secs(jsec)%name)
657            WRITE(numout,*)       "      Compute temperature and salinity transports ? ",secs(jsec)%llstrpond
658            WRITE(numout,*)       "      Compute ice transport ?           ",secs(jsec)%ll_ice_section
659            WRITE(numout,*)       "      Section crosses date-line ?       ",secs(jsec)%ll_date_line
660            WRITE(numout,*)       "      Slope section :                   ",secs(jsec)%slopeSection
661            WRITE(numout,*)       "      Number of points in the section:  ",iptglo
662            WRITE(numout,*)       "      Number of classes                 ",secs(jsec)%nb_class
663            WRITE(numout,clformat)"      Insitu density classes :          ",secs(jsec)%zsigi
664            WRITE(numout,clformat)"      Potential density classes :       ",secs(jsec)%zsigp
665            WRITE(numout,clformat)"      Salinity classes :                ",secs(jsec)%zsal
666            WRITE(numout,clformat)"      Temperature classes :             ",secs(jsec)%ztem
667            WRITE(numout,clformat)"      Depth classes :                   ",secs(jsec)%zlay
668        ENDIF               
669
670        IF( iptglo .NE. 0 )THEN
671             
672           !read points'coordinates and directions
673           !--------------------------------------
674           IF ( ln_NOOS .AND. verbose ) THEN
675              WRITE(numout,*) 'Coords and direction read'
676           ENDIF
677
678           coordtemp(:) = POINT_SECTION(0,0) !list of points read
679           directemp(:) = 0                  !value of directions of each points
680           DO jpt=1,iptglo
681              READ(numdct_in)i1,i2
682              coordtemp(jpt)%I = i1 
683              coordtemp(jpt)%J = i2
684           ENDDO
685           READ(numdct_in)directemp(1:iptglo)
686   
687           !debug
688           !-----
689           IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
690              WRITE(numout,*)"      List of points in global domain:"
691              DO jpt=1,iptglo
692                 WRITE(numout,*)'        # I J ',jpt,coordtemp(jpt),directemp(jpt)
693              ENDDO                 
694           ENDIF
695
696           !Now each proc selects only points that are in its domain:
697           !--------------------------------------------------------
698           iptloc = 0                    !initialize number of points selected
699           DO jpt=1,iptglo               !loop on listpoint read in the file
700                   
701              iiglo=coordtemp(jpt)%I          ! global coordinates of the point
702              ijglo=coordtemp(jpt)%J          !  "
703
704              IF( iiglo==jpidta .AND. nimpp==1 ) iiglo = 2
705
706              iiloc=iiglo-jpizoom+1-nimpp+1   ! local coordinates of the point
707              ijloc=ijglo-jpjzoom+1-njmpp+1   !  "
708
709              !verify if the point is on the local domain:(1,nlei)*(1,nlej)
710              IF( iiloc .GE. 1 .AND. iiloc .LE. nlei .AND. &
711                  ijloc .GE. 1 .AND. ijloc .LE. nlej       )THEN
712                 iptloc = iptloc + 1                                                 ! count local points
713                 secs(jsec)%listPoint(iptloc) = POINT_SECTION(mi0(iiglo),mj0(ijglo)) ! store local coordinates
714                 secs(jsec)%direction(iptloc) = directemp(jpt)                       ! store local direction
715              ENDIF
716
717           ENDDO
718     
719           secs(jsec)%nb_point=iptloc !store number of section's points
720
721           !debug
722           !-----
723           IF(   lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
724              WRITE(numout,*)"      List of points selected by the proc:"
725              DO jpt = 1,iptloc
726                 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
727                 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1
728                 WRITE(numout,*)'         # I J : ',iiglo,ijglo
729              ENDDO
730           ENDIF
731
732           IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN
733              DO jpt = 1,iptloc
734                 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
735                 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1
736              ENDDO
737           ENDIF
738
739           !remove redundant points between processors
740           !------------------------------------------
741           lldebug = .FALSE. ; IF ( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND. lwp ) lldebug = .TRUE.
742           IF( iptloc .NE. 0 )THEN
743              CALL removepoints(secs(jsec),'I','top_list',lldebug)
744              CALL removepoints(secs(jsec),'I','bot_list',lldebug)
745              CALL removepoints(secs(jsec),'J','top_list',lldebug)
746              CALL removepoints(secs(jsec),'J','bot_list',lldebug)
747           ENDIF
748           IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN
749              DO jpt = 1,secs(jsec)%nb_point
750                 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
751                 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1
752              ENDDO
753           ENDIF
754
755           !debug
756           !-----
757           IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )THEN
758              WRITE(numout,*)"      List of points after removepoints:"
759              iptloc = secs(jsec)%nb_point
760              DO jpt = 1,iptloc
761                 iiglo = secs(jsec)%listPoint(jpt)%I + jpizoom - 1 + nimpp - 1
762                 ijglo = secs(jsec)%listPoint(jpt)%J + jpjzoom - 1 + njmpp - 1
763                 WRITE(numout,*)'         # I J : ',iiglo,ijglo
764                 CALL FLUSH(numout)
765              ENDDO
766           ENDIF
767
768        ELSE  ! iptglo = 0
769           IF( lwp .AND. ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) )&
770              WRITE(numout,*)'   No points for this section.'
771        ENDIF
772
773     ENDDO !end of the loop on jsec
774 
775     nb_sec = jsec-1   !number of section read in the file
776
777     IF( lwp .AND. verbose )  WRITE(numout,*)'diadct: read sections: Finished readsec.'
778
779     CALL wrk_dealloc( nb_point_max, directemp )
780     !
781  END SUBROUTINE readsec
782
783  SUBROUTINE removepoints(sec,cdind,cdextr,ld_debug)
784     !!---------------------------------------------------------------------------
785     !!             *** function removepoints
786     !!
787     !!   ** Purpose :: Remove points which are common to 2 procs
788     !!
789     !----------------------------------------------------------------------------
790     !! * arguments
791     TYPE(SECTION),INTENT(INOUT) :: sec
792     CHARACTER(len=1),INTENT(IN) :: cdind   ! = 'I'/'J'
793     CHARACTER(len=8),INTENT(IN) :: cdextr  ! = 'top_list'/'bot_list'
794     LOGICAL,INTENT(IN)          :: ld_debug                     
795
796     !! * Local variables
797     INTEGER :: iextr         ,& !extremity of listpoint that we verify
798                iind          ,& !coord     of listpoint that we verify
799                itest         ,& !indice value of the side of the domain
800                                 !where points could be redundant
801                isgn          ,& ! isgn= 1 : scan listpoint from start to end
802                                 ! isgn=-1 : scan listpoint from end to start
803                istart,iend      !first and last points selected in listpoint
804     INTEGER :: jpoint           !loop on list points
805     INTEGER, POINTER, DIMENSION(:)   :: idirec !contains temporary sec%direction
806     INTEGER, POINTER, DIMENSION(:,:) :: icoord !contains temporary sec%listpoint
807     !----------------------------------------------------------------------------
808     CALL wrk_alloc(    nb_point_max, idirec )
809     CALL wrk_alloc( 2, nb_point_max, icoord )
810
811     IF( ld_debug )WRITE(numout,*)'      -------------------------'
812     IF( ld_debug )WRITE(numout,*)'      removepoints in listpoint'
813
814     !iextr=extremity of list_point that we verify
815     IF      ( cdextr=='bot_list' )THEN ; iextr=1            ; isgn=1
816     ELSE IF ( cdextr=='top_list' )THEN ; iextr=sec%nb_point ; isgn=-1
817     ELSE    ; CALL ctl_stop("removepoints :Wrong value for cdextr")
818     ENDIF
819 
820     !which coordinate shall we verify ?
821     IF      ( cdind=='I' )THEN   ; itest=nlei ; iind=1
822     ELSE IF ( cdind=='J' )THEN   ; itest=nlej ; iind=2
823     ELSE    ; CALL ctl_stop("removepoints :Wrong value for cdind") 
824     ENDIF
825
826     IF( ld_debug )THEN
827        WRITE(numout,*)'      case: coord/list extr/domain side'
828        WRITE(numout,*)'      ', cdind,' ',cdextr,' ',itest
829        WRITE(numout,*)'      Actual number of points: ',sec%nb_point
830     ENDIF
831
832     icoord(1,1:nb_point_max) = sec%listPoint%I
833     icoord(2,1:nb_point_max) = sec%listPoint%J
834     idirec                   = sec%direction
835     sec%listPoint            = POINT_SECTION(0,0)
836     sec%direction            = 0
837
838     jpoint=iextr+isgn
839     DO WHILE( jpoint .GE. 1 .AND. jpoint .LE. sec%nb_point )
840         IF( icoord( iind,jpoint-isgn ) == itest .AND. icoord( iind,jpoint ) == itest )THEN ; jpoint=jpoint+isgn
841         ELSE                                                                               ; EXIT
842         ENDIF
843     ENDDO 
844
845     IF( cdextr=='bot_list')THEN ; istart=jpoint-1 ; iend=sec%nb_point
846     ELSE                        ; istart=1        ; iend=jpoint+1
847     ENDIF
848
849     sec%listPoint(1:1+iend-istart)%I = icoord(1,istart:iend)
850     sec%listPoint(1:1+iend-istart)%J = icoord(2,istart:iend)
851     sec%direction(1:1+iend-istart)   = idirec(istart:iend)
852     sec%nb_point                     = iend-istart+1
853     
854     IF( ld_debug )THEN
855        WRITE(numout,*)'      Number of points after removepoints :',sec%nb_point
856        WRITE(numout,*)'      sec%direction after removepoints :',sec%direction(1:sec%nb_point)
857     ENDIF
858
859     CALL wrk_dealloc(    nb_point_max, idirec )
860     CALL wrk_dealloc( 2, nb_point_max, icoord )
861  END SUBROUTINE removepoints
862
863  SUBROUTINE transport(sec,ld_debug,jsec)
864     !!-------------------------------------------------------------------------------------------
865     !!                     ***  ROUTINE transport  ***
866     !!
867     !!  Purpose ::  Compute the transport for each point in a section
868     !!
869     !!  Method  ::  Loop over each segment, and each vertical level and add the transport
870     !!              Be aware :           
871     !!              One section is a sum of segments
872     !!              One segment is defined by 2 consecutive points in sec%listPoint
873     !!              All points of sec%listPoint are positioned on the F-point of the cell
874     !!
875     !!              There are two loops:                 
876     !!              loop on the segment between 2 nodes
877     !!              loop on the level jk !!
878     !!
879     !! ** Output: Arrays containing the volume,density,salinity,temperature etc
880     !!            transports for each point in a section, summed over each nn_dct.
881     !!
882     !!-------------------------------------------------------------------------------------------
883     !! * Arguments
884     TYPE(SECTION),INTENT(INOUT) :: sec
885     LOGICAL      ,INTENT(IN)    :: ld_debug
886     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section
887   
888     !! * Local variables
889     INTEGER             :: jk,jseg,jclass,    &!loop on level/segment/classes
890                            isgnu  , isgnv      !
891     REAL(wp):: zumid        , zvmid        ,  &!U/V velocity on a cell segment
892                zumid_ice    , zvmid_ice    ,  &!U/V ice velocity
893                zTnorm                          !transport of velocity through one cell's sides
894     REAL(wp):: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point
895
896     TYPE(POINT_SECTION) :: k
897     !!--------------------------------------------------------
898
899     IF( ld_debug )WRITE(numout,*)'      Compute transport (jsec,sec%nb_point,sec%slopeSection) : ', jsec,sec%nb_point,sec%slopeSection
900
901     !---------------------------!
902     !  COMPUTE TRANSPORT        !
903     !---------------------------!
904     IF(sec%nb_point .NE. 0)THEN   
905
906        !----------------------------------------------------------------------------------------------------
907        !----------------------------------------------------------------------------------------------------
908        !----------------------------------------------------------------------------------------------------
909        !
910        !
911        ! ! ! ! JT 1/09/2018 - changing convention. Always direction + is toward left hand of section
912        !
913        !    Making sign of the velocities used to calculate the volume transport a function of direction, not slopesection
914        !    (isgnu, isgnv)
915        !   
916        !    They vary for each segment of the section.
917        !
918        !----------------------------------------------------------------------------------------------------
919        !----------------------------------------------------------------------------------------------------
920        !----------------------------------------------------------------------------------------------------
921        !Compute sign for velocities:
922        !
923        !convention:
924        !   non horizontal section: direction + is toward left hand of section
925        !       horizontal section: direction + is toward north of section
926        !
927        !
928        !       slopeSection < 0     slopeSection > 0       slopeSection=inf            slopeSection=0
929        !       ----------------      -----------------     ---------------             --------------
930        !
931        !   isgnv=1         direction +     
932        !  ______         _____             ______                                                   
933        !        |           //|            |                  |                         direction +   
934        !        | isgnu=1  // |            |isgnu=1           |isgnu=1                     /|\
935        !        |_______  //         ______|    \\            | ---\                        |
936        !               |             | isgnv=-1  \\ |         | ---/ direction +       ____________
937        !               |             |          __\\|         |                   
938        !               |             |     direction +        |                      isgnv=1                                 
939        !                                                     
940        !----------------------------------------------------------------------------------------------------
941
942        IF( ld_debug )write(numout,*)"sec%slopeSection isgnu isgnv ",sec%slopeSection,isgnu,isgnv
943
944        !--------------------------------------!
945        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  !
946        !--------------------------------------!
947        DO jseg=1,MAX(sec%nb_point-1,0)
948           
949           
950           !Compute sign for velocities:
951           
952           !isgnu =  1
953           !isgnv =  1
954           !
955           !changing sign of u and v is dependent on the direction of the section.
956           !isgnu =  1
957           !isgnv =  1
958           !SELECT CASE( sec%direction(jseg) )
959           !CASE(0)  ;   isgnv = -1
960           !CASE(3)  ;   isgnu = -1
961           !END SELECT
962           
963           
964           SELECT CASE( sec%direction(jseg) )
965           CASE(0) 
966              isgnu =  1
967              isgnv = -1
968           CASE(1)
969              isgnu =  1
970              isgnv =  1
971           CASE(2) 
972              isgnu =  1
973              isgnv =  1
974           CASE(3) 
975              isgnu = -1
976              isgnv =  1
977           END SELECT
978           
979           !-------------------------------------------------------------------------------------------
980           ! Select the appropriate coordinate for computing the velocity of the segment
981           ! Corrected by JT 01/09/2018 (#)
982           !
983           !                      CASE(0)                                    Case (2)
984           !                      -------                                    --------
985           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)     
986           !      F(i,j)---------#V(i,j)-------F(i+1,j)                                 |
987           !                     -------->                                              |
988           !                                                                        |   |
989           !                                                                        |   |
990           !                      Case (3)                                          | U(i,j)
991           !                      --------                                          |   |
992           !                                                                        V   |
993           !  listPoint(jseg+1) F(i,j+1)                                                |
994           !                        |                                                   |
995           !                        |                                                   |
996           !                        |                                 listPoint(jseg+1) F(i,j-1)
997           !                   ^    |                                           
998           !                   |    |                                           
999           !                   | U(i,j+1)                                           
1000           !                   |    |                                       Case(1)     
1001           !                   |    |                                       ------     
1002           !                        |                                           
1003           !                        |                 listPoint(jseg+1)             listPoint(jseg)                           
1004           !                        |                 F(i-1,j)----------#V(i-1,j) ------#f(i,j)                           
1005           ! listPoint(jseg)     F(i,j)                                 <-------
1006           !
1007           !-------------------------------------------------------------------------------------------
1008
1009           SELECT CASE( sec%direction(jseg) )
1010           CASE(0)  ;   k = sec%listPoint(jseg)
1011           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
1012           CASE(2)  ;   k = sec%listPoint(jseg)
1013           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
1014           END SELECT
1015
1016           !---------------------------|
1017           !     LOOP ON THE LEVEL     |
1018           !---------------------------|
1019           !Sum of the transport on the vertical
1020           DO jk=1,mbathy(k%I,k%J)
1021
1022              ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
1023              SELECT CASE( sec%direction(jseg) )
1024              CASE(0,1)
1025                 ztn   = interp(k%I,k%J,jk,'V',0 )
1026                 zsn   = interp(k%I,k%J,jk,'V',1 )
1027                 zrhop = interp(k%I,k%J,jk,'V',2)
1028                 zrhoi = interp(k%I,k%J,jk,'V',3)
1029                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1)
1030              CASE(2,3)
1031                 ztn   = interp(k%I,k%J,jk,'U',0)
1032                 zsn   = interp(k%I,k%J,jk,'U',1)
1033                 zrhop = interp(k%I,k%J,jk,'U',2)
1034                 zrhoi = interp(k%I,k%J,jk,'U',3)
1035                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1) 
1036              END SELECT
1037
1038              zfsdep= fsdept(k%I,k%J,jk)
1039 
1040              !compute velocity with the correct direction
1041              SELECT CASE( sec%direction(jseg) )
1042              CASE(0,1) 
1043                 zumid=0.
1044                 zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk)
1045              CASE(2,3)
1046                 zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk)
1047                 zvmid=0.
1048              END SELECT
1049
1050              !zTnorm=transport through one cell;
1051              !velocity* cell's length * cell's thickness
1052              zTnorm=zumid*e2u(k%I,k%J)*  fse3u(k%I,k%J,jk)+     &
1053                     zvmid*e1v(k%I,k%J)*  fse3v(k%I,k%J,jk)
1054
1055#if ! defined key_vvl
1056              !add transport due to free surface
1057              IF( jk==1 )THEN
1058                 zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + &
1059                                   zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk)
1060              ENDIF
1061#endif
1062              !COMPUTE TRANSPORT 
1063
1064              transports_3d(1,jsec,jseg,jk) = transports_3d(1,jsec,jseg,jk) + zTnorm
1065 
1066              IF ( sec%llstrpond ) THEN
1067                 transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk)  + zTnorm * zrhoi
1068                 transports_3d(3,jsec,jseg,jk) = transports_3d(3,jsec,jseg,jk)  + zTnorm * zrhop
1069                 transports_3d(4,jsec,jseg,jk) = transports_3d(4,jsec,jseg,jk)  + zTnorm * 4.e+3_wp * (ztn+273.15) * 1026._wp
1070                 transports_3d(5,jsec,jseg,jk) = transports_3d(5,jsec,jseg,jk)  + zTnorm * 0.001 * zsn * 1026._wp
1071              ENDIF
1072   
1073           ENDDO !end of loop on the level
1074
1075#if defined key_lim2 || defined key_lim3
1076
1077           !ICE CASE   
1078           !------------
1079           IF( sec%ll_ice_section )THEN
1080              SELECT CASE (sec%direction(jseg))
1081              CASE(0)
1082                 zumid_ice = 0
1083                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
1084              CASE(1)
1085                 zumid_ice = 0
1086                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
1087              CASE(2)
1088                 zvmid_ice = 0
1089                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
1090              CASE(3)
1091                 zvmid_ice = 0
1092                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
1093              END SELECT
1094   
1095              zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J)
1096   
1097              transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)*   &
1098                                   (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &
1099                                  *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) +  &
1100                                    hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) &
1101                                   +zice_vol_pos
1102              transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)*   &
1103                                    (1.0 -  frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &
1104                                   +zice_surf_pos
1105   
1106           ENDIF !end of ice case
1107#endif
1108 
1109        ENDDO !end of loop on the segment
1110
1111     ENDIF !end of sec%nb_point =0 case
1112     !
1113  END SUBROUTINE transport
1114 
1115  SUBROUTINE transport_h(sec,ld_debug,jsec)
1116     !!-------------------------------------------------------------------------------------------
1117     !!                     ***  ROUTINE hourly_transport  ***
1118     !!
1119     !!              Exactly as routine transport but for data summed at
1120     !!              each time step and averaged each hour
1121     !!
1122     !!  Purpose ::  Compute the transport for each point in a section
1123     !!
1124     !!  Method  ::  Loop over each segment, and each vertical level and add the transport
1125     !!              Be aware :         
1126     !!              One section is a sum of segments
1127     !!              One segment is defined by 2 consecutive points in sec%listPoint
1128     !!              All points of sec%listPoint are positioned on the F-point of the cell
1129     !!
1130     !!              There are two loops:                 
1131     !!              loop on the segment between 2 nodes
1132     !!              loop on the level jk
1133     !!
1134     !! ** Output: Arrays containing the volume,density,salinity,temperature etc
1135     !!            transports for each point in a section, summed over each nn_dct.
1136     !!
1137     !!-------------------------------------------------------------------------------------------
1138     !! * Arguments
1139     TYPE(SECTION),INTENT(INOUT) :: sec
1140     LOGICAL      ,INTENT(IN)    :: ld_debug
1141     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section
1142   
1143     !! * Local variables
1144     INTEGER             :: jk,jseg,jclass,    &!loop on level/segment/classes
1145                            isgnu  , isgnv      !
1146     REAL(wp):: zumid        , zvmid        ,  &!U/V velocity on a cell segment
1147                zumid_ice    , zvmid_ice    ,  &!U/V ice velocity
1148                zTnorm                          !transport of velocity through one cell's sides
1149     REAL(wp):: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point
1150
1151     TYPE(POINT_SECTION) :: k
1152     !!--------------------------------------------------------
1153
1154     !!NIALL IF( ld_debug )WRITE(numout,*)'      Compute transport'
1155
1156     !---------------------------!
1157     !  COMPUTE TRANSPORT        !
1158     !---------------------------!
1159     IF(sec%nb_point .NE. 0)THEN   
1160
1161        !----------------------------------------------------------------------------------------------------
1162        !----------------------------------------------------------------------------------------------------
1163        !----------------------------------------------------------------------------------------------------
1164        !
1165        !
1166        ! ! ! ! JT 1/09/2018 - changing convention. Always direction + is toward left hand of section
1167        !
1168        !    Making sign of the velocities used to calculate the volume transport a function of direction, not slopesection
1169        !    (isgnu, isgnv)
1170        !   
1171        !    They vary for each segment of the section.
1172        !
1173        !----------------------------------------------------------------------------------------------------
1174        !----------------------------------------------------------------------------------------------------
1175        !----------------------------------------------------------------------------------------------------
1176        !Compute sign for velocities:
1177        !
1178        !convention:
1179        !   non horizontal section: direction + is toward left hand of section
1180        !       horizontal section: direction + is toward north of section
1181        !
1182        !
1183        !       slopeSection < 0     slopeSection > 0       slopeSection=inf            slopeSection=0
1184        !       ----------------      -----------------     ---------------             --------------
1185        !
1186        !   isgnv=1         direction +     
1187        !  ______         _____             ______                                                   
1188        !        |           //|            |                  |                         direction +   
1189        !        | isgnu=1  // |            |isgnu=1           |isgnu=1                     /|\
1190        !        |_______  //         ______|    \\            | ---\                        |
1191        !               |             | isgnv=-1  \\ |         | ---/ direction +       ____________
1192        !               |             |          __\\|         |                   
1193        !               |             |     direction +        |                      isgnv=1                                 
1194        !                                                     
1195        !----------------------------------------------------------------------------------------------------
1196
1197        IF( ld_debug )write(numout,*)"isgnu isgnv ",isgnu,isgnv
1198
1199        !--------------------------------------!
1200        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  !
1201        !--------------------------------------!
1202        DO jseg=1,MAX(sec%nb_point-1,0)
1203           
1204           
1205           !Compute sign for velocities:
1206           
1207           !isgnu =  1
1208           !isgnv =  1
1209           !
1210           ! changing sign of u and v is dependent on the direction of the section.
1211           !isgnu =  1
1212           !isgnv =  1
1213           !SELECT CASE( sec%direction(jseg) )
1214           !CASE(0)  ;   isgnv = -1
1215           !CASE(3)  ;   isgnu = -1
1216           !END SELECT
1217           
1218           
1219           SELECT CASE( sec%direction(jseg) )
1220           CASE(0) 
1221              isgnu =  1
1222              isgnv = -1
1223           CASE(1)
1224              isgnu =  1
1225              isgnv =  1
1226           CASE(2) 
1227              isgnu =  1
1228              isgnv =  1
1229           CASE(3) 
1230              isgnu = -1
1231              isgnv =  1
1232           END SELECT
1233           
1234           !-------------------------------------------------------------------------------------------
1235           ! Select the appropriate coordinate for computing the velocity of the segment
1236           ! Corrected by JT 01/09/2018 (#)
1237           !
1238           !                      CASE(0)                                    Case (2)
1239           !                      -------                                    --------
1240           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)     
1241           !      F(i,j)---------#V(i,j)-------F(i+1,j)                                 |
1242           !                     -------->                                              |
1243           !                                                                        |   |
1244           !                                                                        |   |
1245           !                      Case (3)                                          | U(i,j)
1246           !                      --------                                          |   |
1247           !                                                                        V   |
1248           !  listPoint(jseg+1) F(i,j+1)                                                |
1249           !                        |                                                   |
1250           !                        |                                                   |
1251           !                        |                                 listPoint(jseg+1) F(i,j-1)
1252           !                   ^    |                                           
1253           !                   |    |                                           
1254           !                   | U(i,j+1)                                           
1255           !                   |    |                                       Case(1)     
1256           !                   |    |                                       ------     
1257           !                        |                                           
1258           !                        |                 listPoint(jseg+1)             listPoint(jseg)                           
1259           !                        |                 F(i-1,j)----------#V(i-1,j) ------#f(i,j)                           
1260           ! listPoint(jseg)     F(i,j)                                 <-------
1261           !
1262           !-------------------------------------------------------------------------------------------
1263
1264           SELECT CASE( sec%direction(jseg) )
1265           CASE(0)  ;   k = sec%listPoint(jseg)
1266           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
1267           CASE(2)  ;   k = sec%listPoint(jseg)
1268           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
1269           END SELECT
1270
1271           !---------------------------|
1272           !     LOOP ON THE LEVEL     |
1273           !---------------------------|
1274           !Sum of the transport on the vertical
1275           DO jk=1,mbathy(k%I,k%J)
1276
1277              ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
1278              SELECT CASE( sec%direction(jseg) )
1279              CASE(0,1)
1280                 ztn   = interp(k%I,k%J,jk,'V',0)
1281                 zsn   = interp(k%I,k%J,jk,'V',1)
1282                 zrhop = interp(k%I,k%J,jk,'V',2)
1283                 zrhoi = interp(k%I,k%J,jk,'V',3)
1284                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1)
1285              CASE(2,3)
1286                 ztn   = interp(k%I,k%J,jk,'U',0)
1287                 zsn   = interp(k%I,k%J,jk,'U',1)
1288                 zrhop = interp(k%I,k%J,jk,'U',2)
1289                 zrhoi = interp(k%I,k%J,jk,'U',3)
1290                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1) 
1291              END SELECT
1292
1293              zfsdep= fsdept(k%I,k%J,jk)
1294 
1295              !compute velocity with the correct direction
1296              SELECT CASE( sec%direction(jseg) )
1297              CASE(0,1) 
1298                 zumid=0.
1299                 zvmid=isgnv*vn(k%I,k%J,jk)*vmask(k%I,k%J,jk)
1300              CASE(2,3)
1301                 zumid=isgnu*un(k%I,k%J,jk)*umask(k%I,k%J,jk)
1302                 zvmid=0.
1303              END SELECT
1304
1305              !zTnorm=transport through one cell;
1306              !velocity* cell's length * cell's thickness
1307              zTnorm=zumid*e2u(k%I,k%J)*  fse3u(k%I,k%J,jk)+     &
1308                     zvmid*e1v(k%I,k%J)*  fse3v(k%I,k%J,jk)
1309
1310#if ! defined key_vvl
1311              !add transport due to free surface
1312              IF( jk==1 )THEN
1313                 zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) + &
1314                                   zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk)
1315              ENDIF
1316#endif
1317              !COMPUTE TRANSPORT
1318
1319              transports_3d_h(1,jsec,jseg,jk) = transports_3d_h(1,jsec,jseg,jk) + zTnorm
1320 
1321              IF ( sec%llstrpond ) THEN
1322                 transports_3d_h(2,jsec,jseg,jk) = transports_3d_h(2,jsec,jseg,jk)  + zTnorm * zrhoi
1323                 transports_3d_h(3,jsec,jseg,jk) = transports_3d_h(3,jsec,jseg,jk)  + zTnorm * zrhop
1324                 transports_3d_h(4,jsec,jseg,jk) = transports_3d_h(4,jsec,jseg,jk)  + zTnorm * 4.e+3_wp * (ztn+273.15) * 1026._wp
1325                 transports_3d_h(5,jsec,jseg,jk) = transports_3d_h(5,jsec,jseg,jk)  + zTnorm * 0.001 * zsn * 1026._wp
1326              ENDIF
1327
1328           ENDDO !end of loop on the level
1329
1330#if defined key_lim2 || defined key_lim3
1331
1332           !ICE CASE   
1333           !------------
1334           IF( sec%ll_ice_section )THEN
1335              SELECT CASE (sec%direction(jseg))
1336              CASE(0)
1337                 zumid_ice = 0
1338                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
1339              CASE(1)
1340                 zumid_ice = 0
1341                 zvmid_ice =  isgnv*0.5*(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
1342              CASE(2)
1343                 zvmid_ice = 0
1344                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
1345              CASE(3)
1346                 zvmid_ice = 0
1347                 zumid_ice =  isgnu*0.5*(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
1348              END SELECT
1349   
1350              zTnorm=zumid_ice*e2u(k%I,k%J)+zvmid_ice*e1v(k%I,k%J)
1351   
1352              transports_2d_h(1,jsec,jseg) = transports_2d_h(1,jsec,jseg) + (zTnorm)*   &
1353                                   (1.0 - frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &
1354                                  *(hsnif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J) +  &
1355                                    hicif(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J)) &
1356                                   +zice_vol_pos
1357              transports_2d_h(2,jsec,jseg) = transports_2d_h(2,jsec,jseg) + (zTnorm)*   &
1358                                    (1.0 -  frld(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J))  &
1359                                   +zice_surf_pos
1360   
1361           ENDIF !end of ice case
1362#endif
1363 
1364        ENDDO !end of loop on the segment
1365
1366     ENDIF   !end of sec%nb_point =0 case
1367     !
1368  END SUBROUTINE transport_h
1369 
1370  SUBROUTINE dia_dct_sum(sec,jsec) 
1371     !!-------------------------------------------------------------
1372     !! Purpose: Average the transport over nn_dctwri time steps 
1373     !! and sum over the density/salinity/temperature/depth classes
1374     !!
1375     !! Method:   Sum over relevant grid cells to obtain values 
1376     !!           for each class
1377     !!              There are several loops:                 
1378     !!              loop on the segment between 2 nodes
1379     !!              loop on the level jk
1380     !!              loop on the density/temperature/salinity/level classes
1381     !!              test on the density/temperature/salinity/level
1382     !!
1383     !!  Note:    Transport through a given section is equal to the sum of transports
1384     !!           computed on each proc.
1385     !!           On each proc,transport is equal to the sum of transport computed through
1386     !!           segments linking each point of sec%listPoint  with the next one.   
1387     !!
1388     !!-------------------------------------------------------------
1389     !! * arguments
1390     TYPE(SECTION),INTENT(INOUT) :: sec 
1391     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section
1392 
1393     TYPE(POINT_SECTION) :: k 
1394     INTEGER  :: jk,jseg,jclass                        ! dummy variables for looping on level/segment/classes 
1395     REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point
1396     !!-------------------------------------------------------------
1397
1398
1399     !! Sum the relevant segments to obtain values for each class
1400     IF(sec%nb_point .NE. 0)THEN   
1401
1402        !--------------------------------------!
1403        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  !
1404        !--------------------------------------!
1405        DO jseg=1,MAX(sec%nb_point-1,0)
1406           
1407           !-------------------------------------------------------------------------------------------
1408           ! Select the appropriate coordinate for computing the velocity of the segment
1409           !
1410           !                      CASE(0)                                    Case (2)
1411           !                      -------                                    --------
1412           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)     
1413           !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               |
1414           !                                                                            |
1415           !                                                                            |
1416           !                                                                            |
1417           !                      Case (3)                                            U(i,j)
1418           !                      --------                                              |
1419           !                                                                            |
1420           !  listPoint(jseg+1) F(i,j+1)                                                |
1421           !                        |                                                   |
1422           !                        |                                                   |
1423           !                        |                                 listPoint(jseg+1) F(i,j-1)
1424           !                        |                                           
1425           !                        |                                           
1426           !                     U(i,j+1)                                           
1427           !                        |                                       Case(1)     
1428           !                        |                                       ------     
1429           !                        |                                           
1430           !                        |                 listPoint(jseg+1)             listPoint(jseg)                           
1431           !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                           
1432           ! listPoint(jseg)     F(i,j)
1433           !
1434           !-------------------------------------------------------------------------------------------
1435
1436           SELECT CASE( sec%direction(jseg) )
1437           CASE(0)  ;   k = sec%listPoint(jseg)
1438           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
1439           CASE(2)  ;   k = sec%listPoint(jseg)
1440           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
1441           END SELECT
1442
1443           !---------------------------|
1444           !     LOOP ON THE LEVEL     |
1445           !---------------------------|
1446           !Sum of the transport on the vertical
1447           DO jk=1,mbathy(k%I,k%J)
1448
1449              ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
1450              SELECT CASE( sec%direction(jseg) )
1451              CASE(0,1)
1452                 ztn   = interp(k%I,k%J,jk,'V',0)
1453                 zsn   = interp(k%I,k%J,jk,'V',1)
1454                 zrhop = interp(k%I,k%J,jk,'V',2)
1455                 zrhoi = interp(k%I,k%J,jk,'V',3)
1456                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1)
1457              CASE(2,3)
1458                 ztn   = interp(k%I,k%J,jk,'U',0)
1459                 zsn   = interp(k%I,k%J,jk,'U',1)
1460                 zrhop = interp(k%I,k%J,jk,'U',2)
1461                 zrhoi = interp(k%I,k%J,jk,'U',3)
1462                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1) 
1463              END SELECT
1464
1465              zfsdep= fsdept(k%I,k%J,jk) 
1466 
1467              !-------------------------------
1468              !  LOOP ON THE DENSITY CLASSES |
1469              !-------------------------------
1470              !The computation is made for each density/temperature/salinity/depth class
1471              DO jclass=1,MAX(1,sec%nb_class-1)
1472 
1473                 !----------------------------------------------!
1474                 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL!
1475                 !----------------------------------------------!
1476 
1477                 IF ( (                                                    &
1478                    ((( zrhop .GE. (sec%zsigp(jclass)+1000.  )) .AND.      &
1479                    (   zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR.     &
1480                    ( sec%zsigp(jclass) .EQ. 99.)) .AND.                   &
1481
1482                    ((( zrhoi .GE. (sec%zsigi(jclass) + 1000.  )) .AND.    &
1483                    (   zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR.     &
1484                    ( sec%zsigi(jclass) .EQ. 99.)) .AND.                   &
1485
1486                    ((( zsn .GT. sec%zsal(jclass)) .AND.                   &
1487                    (   zsn .LE. sec%zsal(jclass+1))) .OR.                 &
1488                    ( sec%zsal(jclass) .EQ. 99.)) .AND.                    &
1489
1490                    ((( ztn .GE. sec%ztem(jclass)) .AND.                   &
1491                    (   ztn .LE. sec%ztem(jclass+1))) .OR.                 &
1492                    ( sec%ztem(jclass) .EQ.99.)) .AND.                     &
1493
1494                    ((( zfsdep .GE. sec%zlay(jclass)) .AND.                &
1495                    (   zfsdep .LE. sec%zlay(jclass+1))) .OR.              &
1496                    ( sec%zlay(jclass) .EQ. 99. ))                         &
1497                                                                   ))   THEN
1498
1499                    !SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS
1500                    !----------------------------------------------------------------------------
1501                    IF (transports_3d(1,jsec,jseg,jk) .GE. 0.0) THEN
1502                       sec%transport(1,jclass) = sec%transport(1,jclass)+transports_3d(1,jsec,jseg,jk)
1503                    ELSE
1504                       sec%transport(2,jclass) = sec%transport(2,jclass)+transports_3d(1,jsec,jseg,jk)
1505                    ENDIF
1506                    IF( sec%llstrpond )THEN
1507
1508                       IF( transports_3d(1,jsec,jseg,jk) .NE. 0._wp ) THEN
1509
1510                          IF (transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) .GE. 0.0 ) THEN
1511                             sec%transport(3,jclass) = sec%transport(3,jclass)+transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk)
1512                          ELSE
1513                             sec%transport(4,jclass) = sec%transport(4,jclass)+transports_3d(2,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk)
1514                          ENDIF
1515
1516                          IF ( transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk) .GE. 0.0 ) THEN
1517                             sec%transport(5,jclass) = sec%transport(5,jclass)+transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk)
1518                          ELSE
1519                             sec%transport(6,jclass) = sec%transport(6,jclass)+transports_3d(3,jsec,jseg,jk)/transports_3d(1,jsec,jseg,jk)
1520                          ENDIF
1521
1522                       ENDIF
1523
1524                       IF ( transports_3d(4,jsec,jseg,jk) .GE. 0.0 ) THEN
1525                          sec%transport(7,jclass) = sec%transport(7,jclass)+transports_3d(4,jsec,jseg,jk)
1526                       ELSE
1527                          sec%transport(8,jclass) = sec%transport(8,jclass)+transports_3d(4,jsec,jseg,jk)
1528                       ENDIF
1529
1530                       IF ( transports_3d(5,jsec,jseg,jk) .GE. 0.0 ) THEN
1531                          sec%transport( 9,jclass) = sec%transport( 9,jclass)+transports_3d(5,jsec,jseg,jk)
1532                       ELSE
1533                          sec%transport(10,jclass) = sec%transport(10,jclass)+transports_3d(5,jsec,jseg,jk)
1534                       ENDIF
1535 
1536                    ELSE
1537                       sec%transport( 3,jclass) = 0._wp 
1538                       sec%transport( 4,jclass) = 0._wp 
1539                       sec%transport( 5,jclass) = 0._wp 
1540                       sec%transport( 6,jclass) = 0._wp 
1541                       sec%transport( 7,jclass) = 0._wp
1542                       sec%transport( 8,jclass) = 0._wp
1543                       sec%transport( 9,jclass) = 0._wp
1544                       sec%transport(10,jclass) = 0._wp
1545                    ENDIF
1546 
1547                 ENDIF ! end of test if point is in class
1548   
1549              ENDDO ! end of loop on the classes
1550 
1551           ENDDO ! loop over jk
1552 
1553#if defined key_lim2 || defined key_lim3 
1554 
1555           !ICE CASE     
1556           IF( sec%ll_ice_section )THEN
1557 
1558              IF ( transports_2d(1,jsec,jseg) .GE. 0.0 ) THEN
1559                 sec%transport(11,1) = sec%transport(11,1)+transports_2d(1,jsec,jseg)
1560              ELSE
1561                 sec%transport(12,1) = sec%transport(12,1)+transports_2d(1,jsec,jseg)
1562              ENDIF
1563 
1564              IF ( transports_2d(3,jsec,jseg) .GE. 0.0 ) THEN
1565                 sec%transport(13,1) = sec%transport(13,1)+transports_2d(2,jsec,jseg)
1566              ELSE
1567                 sec%transport(14,1) = sec%transport(14,1)+transports_2d(2,jsec,jseg)
1568              ENDIF
1569 
1570           ENDIF !end of ice case
1571#endif
1572 
1573        ENDDO !end of loop on the segment
1574
1575     ELSE  !if sec%nb_point =0
1576        sec%transport(1:2,:)=0.
1577        IF (sec%llstrpond) sec%transport(3:10,:)=0.
1578        IF (sec%ll_ice_section) sec%transport( 11:14,:)=0.
1579     ENDIF !end of sec%nb_point =0 case
1580
1581  END SUBROUTINE dia_dct_sum
1582 
1583  SUBROUTINE dia_dct_sum_h(sec,jsec)
1584     !!-------------------------------------------------------------
1585     !! Exactly as dia_dct_sum but for hourly files containing data summed at each time step
1586     !!
1587     !! Purpose: Average the transport over nn_dctwri time steps
1588     !! and sum over the density/salinity/temperature/depth classes
1589     !!
1590     !! Method:
1591     !!           Sum over relevant grid cells to obtain values
1592     !!           for each
1593     !!              There are several loops:                 
1594     !!              loop on the segment between 2 nodes
1595     !!              loop on the level jk
1596     !!              loop on the density/temperature/salinity/level classes
1597     !!              test on the density/temperature/salinity/level
1598     !!
1599     !!  ** Method  :Transport through a given section is equal to the sum of transports
1600     !!              computed on each proc.
1601     !!              On each proc,transport is equal to the sum of transport computed through
1602     !!              segments linking each point of sec%listPoint  with the next one.   
1603     !!
1604     !!-------------------------------------------------------------
1605     !! * arguments
1606     TYPE(SECTION),INTENT(INOUT) :: sec
1607     INTEGER      ,INTENT(IN)    :: jsec        ! numeric identifier of section
1608
1609     TYPE(POINT_SECTION) :: k
1610     INTEGER  :: jk,jseg,jclass                        !loop on level/segment/classes
1611     REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zfsdep ! temperature/salinity/ssh/potential density /depth at u/v point
1612     !!-------------------------------------------------------------
1613
1614     !! Sum the relevant segments to obtain values for each class
1615     IF(sec%nb_point .NE. 0)THEN   
1616
1617        !--------------------------------------!
1618        ! LOOP ON THE SEGMENT BETWEEN 2 NODES  !
1619        !--------------------------------------!
1620        DO jseg=1,MAX(sec%nb_point-1,0)
1621           
1622           !-------------------------------------------------------------------------------------------
1623           ! Select the appropriate coordinate for computing the velocity of the segment
1624           !
1625           !                      CASE(0)                                    Case (2)
1626           !                      -------                                    --------
1627           !  listPoint(jseg)                 listPoint(jseg+1)       listPoint(jseg)  F(i,j)     
1628           !      F(i,j)----------V(i+1,j)-------F(i+1,j)                               |
1629           !                                                                            |
1630           !                                                                            |
1631           !                                                                            |
1632           !                      Case (3)                                            U(i,j)
1633           !                      --------                                              |
1634           !                                                                            |
1635           !  listPoint(jseg+1) F(i,j+1)                                                |
1636           !                        |                                                   |
1637           !                        |                                                   |
1638           !                        |                                 listPoint(jseg+1) F(i,j-1)
1639           !                        |                                           
1640           !                        |                                           
1641           !                     U(i,j+1)                                           
1642           !                        |                                       Case(1)     
1643           !                        |                                       ------     
1644           !                        |                                           
1645           !                        |                 listPoint(jseg+1)             listPoint(jseg)                           
1646           !                        |                 F(i-1,j)-----------V(i,j) -------f(jseg)                           
1647           ! listPoint(jseg)     F(i,j)
1648           !
1649           !-------------------------------------------------------------------------------------------
1650
1651           SELECT CASE( sec%direction(jseg) )
1652           CASE(0)  ;   k = sec%listPoint(jseg)
1653           CASE(1)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
1654           CASE(2)  ;   k = sec%listPoint(jseg)
1655           CASE(3)  ;   k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
1656           END SELECT
1657
1658           !---------------------------|
1659           !     LOOP ON THE LEVEL     |
1660           !---------------------------|
1661           !Sum of the transport on the vertical
1662           DO jk=1,mbathy(k%I,k%J)
1663
1664              ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
1665              SELECT CASE( sec%direction(jseg) )
1666              CASE(0,1)
1667                 ztn   = interp(k%I,k%J,jk,'V',0)
1668                 zsn   = interp(k%I,k%J,jk,'V',1)
1669                 zrhop = interp(k%I,k%J,jk,'V',2)
1670                 zrhoi = interp(k%I,k%J,jk,'V',3)
1671                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I,k%J+1)    ) * vmask(k%I,k%J,1)
1672              CASE(2,3)
1673                 ztn   = interp(k%I,k%J,jk,'U',0)
1674                 zsn   = interp(k%I,k%J,jk,'U',1)
1675                 zrhop = interp(k%I,k%J,jk,'U',2)
1676                 zrhoi = interp(k%I,k%J,jk,'U',3)
1677                 zsshn =  0.5*( sshn(k%I,k%J)    + sshn(k%I+1,k%J)    ) * umask(k%I,k%J,1) 
1678              END SELECT
1679
1680              zfsdep= fsdept(k%I,k%J,jk)
1681 
1682              !-------------------------------
1683              !  LOOP ON THE DENSITY CLASSES |
1684              !-------------------------------
1685              !The computation is made for each density/heat/salt/... class
1686              DO jclass=1,MAX(1,sec%nb_class-1)
1687
1688                 !----------------------------------------------!
1689                 !TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL!
1690                 !----------------------------------------------!
1691 
1692                 IF ( (                                                    &
1693                    ((( zrhop .GE. (sec%zsigp(jclass)+1000.  )) .AND.      &
1694                    (   zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR.     &
1695                    ( sec%zsigp(jclass) .EQ. 99.)) .AND.                   &
1696
1697                    ((( zrhoi .GE. (sec%zsigi(jclass) + 1000.  )) .AND.    &
1698                    (   zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR.     &
1699                    ( sec%zsigi(jclass) .EQ. 99.)) .AND.                   &
1700
1701                    ((( zsn .GT. sec%zsal(jclass)) .AND.                   &
1702                    (   zsn .LE. sec%zsal(jclass+1))) .OR.                 &
1703                    ( sec%zsal(jclass) .EQ. 99.)) .AND.                    &
1704
1705                    ((( ztn .GE. sec%ztem(jclass)) .AND.                   &
1706                    (   ztn .LE. sec%ztem(jclass+1))) .OR.                 &
1707                    ( sec%ztem(jclass) .EQ.99.)) .AND.                     &
1708
1709                    ((( zfsdep .GE. sec%zlay(jclass)) .AND.                &
1710                    (   zfsdep .LE. sec%zlay(jclass+1))) .OR.              &
1711                    ( sec%zlay(jclass) .EQ. 99. ))                         &
1712                                                                   ))   THEN
1713
1714                    !SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS
1715                    !----------------------------------------------------------------------------
1716                    IF (transports_3d_h(1,jsec,jseg,jk) .GE. 0.0) THEN
1717                       sec%transport_h(1,jclass) = sec%transport_h(1,jclass)+transports_3d_h(1,jsec,jseg,jk)
1718                    ELSE
1719                       sec%transport_h(2,jclass) = sec%transport_h(2,jclass)+transports_3d_h(1,jsec,jseg,jk)
1720                    ENDIF
1721                    IF( sec%llstrpond )THEN
1722
1723                       IF( transports_3d_h(1,jsec,jseg,jk) .NE. 0._wp ) THEN
1724
1725                          IF (transports_3d_h(2,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) .GE. 0.0 ) THEN
1726                             sec%transport_h(3,jclass) = sec%transport_h(3,jclass)+transports_3d_h(2,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk)
1727                          ELSE
1728                             sec%transport_h(4,jclass) = sec%transport_h(4,jclass)+transports_3d_h(2,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk)
1729                          ENDIF
1730
1731                          IF ( transports_3d_h(3,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk) .GE. 0.0 ) THEN
1732                             sec%transport_h(5,jclass) = sec%transport_h(5,jclass)+transports_3d_h(3,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk)
1733                          ELSE
1734                             sec%transport_h(6,jclass) = sec%transport_h(6,jclass)+transports_3d_h(3,jsec,jseg,jk)/transports_3d_h(1,jsec,jseg,jk)
1735                          ENDIF
1736
1737                       ENDIF
1738
1739                       IF ( transports_3d_h(4,jsec,jseg,jk) .GE. 0.0 ) THEN
1740                          sec%transport_h(7,jclass) = sec%transport_h(7,jclass)+transports_3d_h(4,jsec,jseg,jk)
1741                       ELSE
1742                          sec%transport_h(8,jclass) = sec%transport_h(8,jclass)+transports_3d_h(4,jsec,jseg,jk)
1743                       ENDIF
1744
1745                       IF ( transports_3d_h(5,jsec,jseg,jk) .GE. 0.0 ) THEN
1746                          sec%transport_h( 9,jclass) = sec%transport_h( 9,jclass)+transports_3d_h(5,jsec,jseg,jk)
1747                       ELSE
1748                          sec%transport_h(10,jclass) = sec%transport_h(10,jclass)+transports_3d_h(5,jsec,jseg,jk)
1749                       ENDIF
1750
1751                    ELSE
1752                       sec%transport_h( 3,jclass) = 0._wp
1753                       sec%transport_h( 4,jclass) = 0._wp
1754                       sec%transport_h( 5,jclass) = 0._wp
1755                       sec%transport_h( 6,jclass) = 0._wp
1756                       sec%transport_h( 7,jclass) = 0._wp
1757                       sec%transport_h( 8,jclass) = 0._wp
1758                       sec%transport_h( 9,jclass) = 0._wp
1759                       sec%transport_h(10,jclass) = 0._wp
1760                    ENDIF
1761
1762                 ENDIF ! end of test if point is in class
1763   
1764              ENDDO ! end of loop on the classes
1765
1766           ENDDO ! loop over jk
1767
1768#if defined key_lim2 || defined key_lim3
1769
1770           !ICE CASE   
1771           IF( sec%ll_ice_section )THEN
1772
1773              IF ( transports_2d_h(1,jsec,jseg) .GE. 0.0 ) THEN
1774                 sec%transport_h(11,1) = sec%transport_h(11,1)+transports_2d_h(1,jsec,jseg)
1775              ELSE
1776                 sec%transport_h(12,1) = sec%transport_h(12,1)+transports_2d_h(1,jsec,jseg)
1777              ENDIF
1778
1779              IF ( transports_2d_h(3,jsec,jseg) .GE. 0.0 ) THEN
1780                 sec%transport_h(13,1) = sec%transport_h(13,1)+transports_2d_h(2,jsec,jseg)
1781              ELSE
1782                 sec%transport_h(14,1) = sec%transport_h(14,1)+transports_2d_h(2,jsec,jseg)
1783              ENDIF
1784
1785           ENDIF !end of ice case
1786#endif
1787 
1788        ENDDO !end of loop on the segment
1789
1790     ELSE  !if sec%nb_point =0
1791        sec%transport_h(1:2,:)=0.
1792        IF (sec%llstrpond) sec%transport_h(3:10,:)=0.
1793        IF (sec%ll_ice_section) sec%transport_h( 11:14,:)=0.
1794     ENDIF !end of sec%nb_point =0 case
1795
1796  END SUBROUTINE dia_dct_sum_h
1797 
1798  SUBROUTINE dia_dct_wri_NOOS(kt,ksec,sec)
1799     !!-------------------------------------------------------------
1800     !! Write transport output in numdct using NOOS formatting
1801     !!
1802     !! Purpose: Write  transports in ascii files
1803     !!
1804     !! Method:
1805     !!        1. Write volume transports in "volume_transport"
1806     !!           Unit: Sv : area * Velocity / 1.e6
1807     !!
1808     !!        2. Write heat transports in "heat_transport"
1809     !!           Unit: Peta W : area * Velocity * T * rhau * Cp / 1.e15
1810     !!
1811     !!        3. Write salt transports in "salt_transport"
1812     !!           Unit: 10^9 g m^3 / s : area * Velocity * S / 1.e6
1813     !!
1814     !!-------------------------------------------------------------
1815     !!arguments
1816     INTEGER, INTENT(IN)          :: kt          ! time-step
1817     TYPE(SECTION), INTENT(INOUT) :: sec         ! section to write   
1818     INTEGER ,INTENT(IN)          :: ksec        ! section number
1819
1820     !!local declarations
1821     INTEGER               :: jclass,ji             ! Dummy loop
1822     CHARACTER(len=2)      :: classe             ! Classname
1823     REAL(wp)              :: zbnd1,zbnd2        ! Class bounds
1824     REAL(wp)              :: zslope             ! section's slope coeff
1825     !
1826     REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace
1827     CHARACTER(len=3)      :: noos_sect_name             ! Classname
1828     CHARACTER(len=25)      :: noos_var_sect_name
1829     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   noos_iom_dummy
1830     INTEGER               :: IERR
1831     
1832     REAL(wp), DIMENSION(3) :: tmp_iom_output
1833     REAL(wp)               :: max_iom_val
1834     LOGICAL       ::   verbose     
1835     verbose = .false.
1836     
1837     !!-------------------------------------------------------------
1838     
1839     
1840     
1841     IF( lwp .AND. verbose ) THEN
1842        WRITE(numout,*) " "
1843        WRITE(numout,*) "dia_dct_wri_NOOS: write transports through sections at timestep: ", kt
1844        WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
1845     ENDIF   
1846       
1847     CALL wrk_alloc(nb_type_class , zsumclasses ) 
1848
1849     zsumclasses(:)=0._wp
1850     zslope = sec%slopeSection       
1851     
1852     IF( lwp ) THEN
1853         IF  ( ln_dct_ascii ) THEN
1854             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
1855         ELSE
1856             WRITE(numdct_NOOS) nyear,nmonth,nday,ksec-1,sec%nb_class-1,sec%name
1857         ENDIF 
1858     ENDIF
1859   
1860     ! Sum all classes together, to give one values per type (pos tran, neg vol trans etc...).
1861     DO jclass=1,MAX(1,sec%nb_class-1)
1862        zsumclasses(1:nb_type_class)=zsumclasses(1:nb_type_class)+sec%transport(1:nb_type_class,jclass)
1863     ENDDO
1864 
1865     classe   = 'total   '
1866     zbnd1   = 0._wp
1867     zbnd2   = 0._wp
1868     
1869     
1870     
1871     write (noos_sect_name, "(I0.3)")  ksec
1872     
1873     IF ( ln_dct_iom_cont ) THEN
1874         max_iom_val = 1.e10
1875         ALLOCATE( noos_iom_dummy(jpi,jpj,3),  STAT= ierr )
1876            IF( ierr /= 0 )   CALL ctl_stop( 'dia_dct_wri_NOOS: failed to allocate noos_iom_dummy array' )
1877     ENDIF
1878     
1879!     JT   I think changing the sign on the output based on the zslope value is redunant.
1880!
1881!     IF ( zslope .gt. 0._wp .and. zslope .ne. 10000._wp ) THEN
1882!         
1883!         IF( lwp ) THEN
1884!             WRITE(numdct_NOOS,'(9e12.4E2)') -(zsumclasses( 1)+zsumclasses( 2)), -zsumclasses( 2),-zsumclasses( 1),   &
1885!                                            -(zsumclasses( 7)+zsumclasses( 8)), -zsumclasses( 8),-zsumclasses( 7),   &
1886!                                            -(zsumclasses( 9)+zsumclasses(10)), -zsumclasses(10),-zsumclasses( 9)
1887!             CALL FLUSH(numdct_NOOS)
1888!         endif
1889
1890!         
1891!         IF ( ln_dct_iom_cont ) THEN
1892!         
1893!             noos_iom_dummy(:,:,:) = 0.
1894!             
1895!             tmp_iom_output(:) = 0.
1896!             tmp_iom_output(1) = -(zsumclasses( 1)+zsumclasses( 2))
1897!             tmp_iom_output(2) = -zsumclasses( 2)
1898!             tmp_iom_output(3) = -zsumclasses( 1)
1899!             
1900!             ! Convert to Sv
1901!             tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
1902!             tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
1903!             tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
1904!             
1905!             ! limit maximum and minimum values in iom_put
1906!             if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
1907!             if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
1908!             if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
1909!             if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
1910!             if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
1911!             if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
1912!             
1913!             ! Set NaN's to Zero         
1914!             if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
1915!             if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
1916!             if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
1917!             
1918!             noos_iom_dummy(:,:,1) = tmp_iom_output(1)
1919!             noos_iom_dummy(:,:,2) = tmp_iom_output(2)
1920!             noos_iom_dummy(:,:,3) = tmp_iom_output(3)
1921!             
1922!             noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_trans'
1923!             if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
1924!             CALL iom_put( noos_var_sect_name,  noos_iom_dummy )
1925!             noos_iom_dummy(:,:,:) = 0.         
1926!             tmp_iom_output(:) = 0.
1927!             tmp_iom_output(1) = -(zsumclasses( 7)+zsumclasses( 8))
1928!             tmp_iom_output(2) = -zsumclasses( 8)
1929!             tmp_iom_output(3) = -zsumclasses( 7)
1930!             
1931!             ! Convert to TJ/s
1932!             tmp_iom_output(1) = tmp_iom_output(1)*1.E-12
1933!             tmp_iom_output(2) = tmp_iom_output(2)*1.E-12
1934!             tmp_iom_output(3) = tmp_iom_output(3)*1.E-12
1935!             
1936!             ! limit maximum and minimum values in iom_put
1937!             if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
1938!             if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
1939!             if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
1940!             if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
1941!             if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
1942!             if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
1943!             
1944!             ! Set NaN's to Zero         
1945!             if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
1946!             if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
1947!             if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
1948!             
1949!             noos_iom_dummy(:,:,1) = tmp_iom_output(1)
1950!             noos_iom_dummy(:,:,2) = tmp_iom_output(2)
1951!             noos_iom_dummy(:,:,3) = tmp_iom_output(3)
1952!             
1953!             !noos_iom_dummy(:,:,1) = -(zsumclasses( 7)+zsumclasses( 8))
1954!             !noos_iom_dummy(:,:,2) = -zsumclasses( 8)
1955!             !noos_iom_dummy(:,:,3) = -zsumclasses( 7)
1956!             
1957!             noos_var_sect_name =  "noos_" // trim(noos_sect_name) // '_heat'
1958!             if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
1959!             CALL iom_put( noos_var_sect_name,  noos_iom_dummy )
1960!             
1961!             noos_iom_dummy(:,:,:) = 0.         
1962!             tmp_iom_output(:) = 0.         
1963!             tmp_iom_output(1) = -(zsumclasses( 9)+zsumclasses( 10))
1964!             tmp_iom_output(2) = -zsumclasses( 10)
1965!             tmp_iom_output(3) = -zsumclasses( 9)
1966!             
1967!             ! Convert to  MT/s
1968!             tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
1969!             tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
1970!             tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
1971!             
1972!             ! limit maximum and minimum values in iom_put
1973!             if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
1974!             if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
1975!             if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
1976!             if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
1977!             if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
1978!             if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
1979!             
1980!             ! Set NaN's to Zero         
1981!             if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
1982!             if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
1983!             if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
1984!             
1985!             noos_iom_dummy(:,:,1) = tmp_iom_output(1)
1986!             noos_iom_dummy(:,:,2) = tmp_iom_output(2)
1987!             noos_iom_dummy(:,:,3) = tmp_iom_output(3)
1988!             
1989!             !noos_iom_dummy(:,:,1) = -(zsumclasses( 9)+zsumclasses( 10))
1990!             !noos_iom_dummy(:,:,2) = -zsumclasses( 10)
1991!             !noos_iom_dummy(:,:,3) = -zsumclasses( 9)
1992!             
1993!             noos_var_sect_name =  "noos_" // trim(noos_sect_name) // '_salt'
1994!             if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
1995!             CALL iom_put( noos_var_sect_name,  noos_iom_dummy )
1996!             noos_iom_dummy(:,:,:) = 0.     
1997!             tmp_iom_output(:) = 0.   
1998!         ENDIF
1999!     ELSE
2000!        IF( lwp ) THEN
2001!            WRITE(numdct_NOOS,'(9e12.4E2)')   zsumclasses( 1)+zsumclasses( 2) ,  zsumclasses( 1), zsumclasses( 2),   &
2002!                                            zsumclasses( 7)+zsumclasses( 8) ,  zsumclasses( 7), zsumclasses( 8),   &
2003!                                            zsumclasses( 9)+zsumclasses(10) ,  zsumclasses( 9), zsumclasses(10)
2004!            CALL FLUSH(numdct_NOOS)
2005!        endif
2006!         
2007!         
2008!        IF ( ln_dct_iom_cont ) THEN
2009!           
2010!             noos_iom_dummy(:,:,:) = 0.
2011!             tmp_iom_output(:) = 0.
2012!             
2013!             tmp_iom_output(1) = (zsumclasses( 1)+zsumclasses( 2))
2014!             tmp_iom_output(2) = zsumclasses( 1)
2015!             tmp_iom_output(3) = zsumclasses( 2)
2016!             
2017!             ! Convert to Sv
2018!             tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
2019!             tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
2020!             tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
2021!             
2022!             ! limit maximum and minimum values in iom_put
2023!             if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
2024!             if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
2025!             if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
2026!             if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
2027!             if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
2028!             if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
2029!             
2030!             ! Set NaN's to Zero         
2031!             if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
2032!             if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
2033!             if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
2034!             
2035!             noos_iom_dummy(:,:,1) = tmp_iom_output(1)
2036!             noos_iom_dummy(:,:,2) = tmp_iom_output(2)
2037!             noos_iom_dummy(:,:,3) = tmp_iom_output(3)
2038!             
2039!             !noos_iom_dummy(:,:,1) = (zsumclasses( 1)+zsumclasses( 2))
2040!             !noos_iom_dummy(:,:,2) = zsumclasses( 1)
2041!             !noos_iom_dummy(:,:,3) = zsumclasses( 2)
2042!             
2043!             
2044!             
2045!             noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_trans'
2046!             if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
2047!             CALL iom_put( noos_var_sect_name,  noos_iom_dummy )
2048!             noos_iom_dummy(:,:,:) = 0.
2049!             tmp_iom_output(:) = 0.
2050!             
2051!             tmp_iom_output(1) = (zsumclasses( 7)+zsumclasses( 8))
2052!             tmp_iom_output(2) = zsumclasses( 7)
2053!             tmp_iom_output(3) = zsumclasses( 8)
2054!             
2055!             ! Convert to TJ/s
2056!             tmp_iom_output(1) = tmp_iom_output(1)*1.E-12
2057!             tmp_iom_output(2) = tmp_iom_output(2)*1.E-12
2058!             tmp_iom_output(3) = tmp_iom_output(3)*1.E-12
2059!             
2060!             ! limit maximum and minimum values in iom_put
2061!             if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
2062!             if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
2063!             if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
2064!             if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
2065!             if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
2066!             if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
2067!             
2068!             ! Set NaN's to Zero         
2069!             if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
2070!             if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
2071!             if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
2072!             
2073!             noos_iom_dummy(:,:,1) = tmp_iom_output(1)
2074!             noos_iom_dummy(:,:,2) = tmp_iom_output(2)
2075!             noos_iom_dummy(:,:,3) = tmp_iom_output(3)
2076!             
2077!             !noos_iom_dummy(:,:,1) = (zsumclasses( 7)+zsumclasses( 8))
2078!             !noos_iom_dummy(:,:,2) = zsumclasses( 7)
2079!             !noos_iom_dummy(:,:,3) = zsumclasses( 8)
2080!             
2081!             noos_var_sect_name =  "noos_" // trim(noos_sect_name) // '_heat'
2082!             if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
2083!             CALL iom_put(noos_var_sect_name,  noos_iom_dummy ) 
2084!             noos_iom_dummy(:,:,:) = 0.
2085!             tmp_iom_output(:) = 0.
2086!             
2087!             tmp_iom_output(1) = (zsumclasses( 9)+zsumclasses( 10))
2088!             tmp_iom_output(2) = zsumclasses( 9)
2089!             tmp_iom_output(3) = zsumclasses( 10)
2090!             
2091!             ! Convert to  MT/s
2092!             tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
2093!             tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
2094!             tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
2095!             
2096!             
2097!             ! limit maximum and minimum values in iom_put
2098!             if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
2099!             if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
2100!             if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
2101!             if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
2102!             if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
2103!             if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
2104!             
2105!             ! Set NaN's to Zero         
2106!             if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
2107!             if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
2108!             if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
2109!             
2110!             noos_iom_dummy(:,:,1) = tmp_iom_output(1)
2111!             noos_iom_dummy(:,:,2) = tmp_iom_output(2)
2112!             noos_iom_dummy(:,:,3) = tmp_iom_output(3)
2113!             
2114!             !noos_iom_dummy(:,:,1) = (zsumclasses( 9)+zsumclasses( 10))
2115!             !noos_iom_dummy(:,:,2) = zsumclasses( 9)
2116!             !noos_iom_dummy(:,:,3) = zsumclasses( 10)
2117!             
2118!             noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_salt'
2119!             if ( lwp ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
2120!             CALL iom_put(noos_var_sect_name,  noos_iom_dummy )
2121!             noos_iom_dummy(:,:,:) = 0.         
2122!             tmp_iom_output(:) = 0.
2123!         ENDIF
2124!         
2125!     ENDIF
2126     
2127
2128
2129
2130
2131
2132
2133
2134     
2135     
2136     IF( lwp ) THEN
2137        IF  ( ln_dct_ascii ) THEN
2138             !WRITE(numdct_NOOS,'(9e12.4E2)')   zsumclasses( 1)+zsumclasses( 2) ,  zsumclasses( 1), zsumclasses( 2),   &
2139             WRITE(numdct_NOOS,'(3F18.3,6e16.8E2)')   zsumclasses( 1)+zsumclasses( 2) ,  zsumclasses( 1), zsumclasses( 2),   &
2140                                             zsumclasses( 7)+zsumclasses( 8) ,  zsumclasses( 7), zsumclasses( 8),   &
2141                                             zsumclasses( 9)+zsumclasses(10) ,  zsumclasses( 9), zsumclasses(10)
2142             CALL FLUSH(numdct_NOOS)
2143        ELSE
2144             WRITE(numdct_NOOS)   zsumclasses( 1)+zsumclasses( 2) ,  zsumclasses( 1), zsumclasses( 2),   &
2145                                  zsumclasses( 7)+zsumclasses( 8) ,  zsumclasses( 7), zsumclasses( 8),   &
2146                                  zsumclasses( 9)+zsumclasses(10) ,  zsumclasses( 9), zsumclasses(10)
2147             CALL FLUSH(numdct_NOOS) 
2148         ENDIF
2149     ENDIF
2150         
2151    IF ( ln_dct_iom_cont ) THEN
2152       
2153         noos_iom_dummy(:,:,:) = 0.
2154         tmp_iom_output(:) = 0.
2155         
2156         tmp_iom_output(1) = (zsumclasses( 1)+zsumclasses( 2))
2157         tmp_iom_output(2) = zsumclasses( 1)
2158         tmp_iom_output(3) = zsumclasses( 2)
2159         
2160         ! Convert to Sv
2161         tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
2162         tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
2163         tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
2164         
2165         ! limit maximum and minimum values in iom_put
2166         if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
2167         if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
2168         if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
2169         if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
2170         if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
2171         if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
2172         
2173         ! Set NaN's to Zero         
2174         if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
2175         if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
2176         if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
2177         
2178         noos_iom_dummy(:,:,1) = tmp_iom_output(1)
2179         noos_iom_dummy(:,:,2) = tmp_iom_output(2)
2180         noos_iom_dummy(:,:,3) = tmp_iom_output(3)
2181         
2182         !noos_iom_dummy(:,:,1) = (zsumclasses( 1)+zsumclasses( 2))
2183         !noos_iom_dummy(:,:,2) = zsumclasses( 1)
2184         !noos_iom_dummy(:,:,3) = zsumclasses( 2)
2185         
2186         
2187         
2188         noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_trans'
2189         if ( lwp .AND. verbose ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
2190         CALL iom_put( noos_var_sect_name,  noos_iom_dummy )
2191         noos_iom_dummy(:,:,:) = 0.
2192         tmp_iom_output(:) = 0.
2193         
2194         tmp_iom_output(1) = (zsumclasses( 7)+zsumclasses( 8))
2195         tmp_iom_output(2) = zsumclasses( 7)
2196         tmp_iom_output(3) = zsumclasses( 8)
2197         
2198         ! Convert to TJ/s
2199         tmp_iom_output(1) = tmp_iom_output(1)*1.E-12
2200         tmp_iom_output(2) = tmp_iom_output(2)*1.E-12
2201         tmp_iom_output(3) = tmp_iom_output(3)*1.E-12
2202         
2203         ! limit maximum and minimum values in iom_put
2204         if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
2205         if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
2206         if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
2207         if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
2208         if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
2209         if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
2210         
2211         ! Set NaN's to Zero         
2212         if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
2213         if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
2214         if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
2215         
2216         noos_iom_dummy(:,:,1) = tmp_iom_output(1)
2217         noos_iom_dummy(:,:,2) = tmp_iom_output(2)
2218         noos_iom_dummy(:,:,3) = tmp_iom_output(3)
2219         
2220         !noos_iom_dummy(:,:,1) = (zsumclasses( 7)+zsumclasses( 8))
2221         !noos_iom_dummy(:,:,2) = zsumclasses( 7)
2222         !noos_iom_dummy(:,:,3) = zsumclasses( 8)
2223         
2224         noos_var_sect_name =  "noos_" // trim(noos_sect_name) // '_heat'
2225         if ( lwp .AND. verbose ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
2226         CALL iom_put(noos_var_sect_name,  noos_iom_dummy ) 
2227         noos_iom_dummy(:,:,:) = 0.
2228         tmp_iom_output(:) = 0.
2229         
2230         tmp_iom_output(1) = (zsumclasses( 9)+zsumclasses( 10))
2231         tmp_iom_output(2) = zsumclasses( 9)
2232         tmp_iom_output(3) = zsumclasses( 10)
2233         
2234         ! Convert to  MT/s
2235         tmp_iom_output(1) = tmp_iom_output(1)*1.E-6
2236         tmp_iom_output(2) = tmp_iom_output(2)*1.E-6
2237         tmp_iom_output(3) = tmp_iom_output(3)*1.E-6
2238         
2239         
2240         ! limit maximum and minimum values in iom_put
2241         if ( tmp_iom_output(1) .gt.  max_iom_val ) tmp_iom_output(1) =  max_iom_val
2242         if ( tmp_iom_output(1) .lt. -max_iom_val ) tmp_iom_output(1) = -max_iom_val
2243         if ( tmp_iom_output(2) .gt.  max_iom_val ) tmp_iom_output(2) =  max_iom_val
2244         if ( tmp_iom_output(2) .lt. -max_iom_val ) tmp_iom_output(2) = -max_iom_val
2245         if ( tmp_iom_output(3) .gt.  max_iom_val ) tmp_iom_output(3) =  max_iom_val
2246         if ( tmp_iom_output(3) .lt. -max_iom_val ) tmp_iom_output(3) = -max_iom_val
2247         
2248         ! Set NaN's to Zero         
2249         if ( tmp_iom_output(1) .ne.  tmp_iom_output(1) ) tmp_iom_output(1) =  max_iom_val*2
2250         if ( tmp_iom_output(2) .ne.  tmp_iom_output(2) ) tmp_iom_output(1) =  max_iom_val*2
2251         if ( tmp_iom_output(3) .ne.  tmp_iom_output(3) ) tmp_iom_output(1) =  max_iom_val*2
2252         
2253         noos_iom_dummy(:,:,1) = tmp_iom_output(1)
2254         noos_iom_dummy(:,:,2) = tmp_iom_output(2)
2255         noos_iom_dummy(:,:,3) = tmp_iom_output(3)
2256         
2257         !noos_iom_dummy(:,:,1) = (zsumclasses( 9)+zsumclasses( 10))
2258         !noos_iom_dummy(:,:,2) = zsumclasses( 9)
2259         !noos_iom_dummy(:,:,3) = zsumclasses( 10)
2260         
2261         noos_var_sect_name = "noos_" // trim(noos_sect_name) // '_salt'
2262         if ( lwp .AND. verbose ) WRITE(numout,*) 'dia_dct_wri_NOOS iom_put: ', kt,ksec, noos_var_sect_name
2263         CALL iom_put(noos_var_sect_name,  noos_iom_dummy )
2264         noos_iom_dummy(:,:,:) = 0.         
2265         tmp_iom_output(:) = 0.
2266
2267
2268        DEALLOCATE(noos_iom_dummy)
2269     ENDIF
2270     
2271
2272     DO jclass=1,MAX(1,sec%nb_class-1)
2273   
2274        classe   = 'N       '
2275        zbnd1   = 0._wp
2276        zbnd2   = 0._wp
2277
2278        !insitu density classes transports
2279        IF( ( sec%zsigi(jclass)   .NE. 99._wp ) .AND. &
2280            ( sec%zsigi(jclass+1) .NE. 99._wp )       )THEN
2281           classe = 'DI       '
2282           zbnd1 = sec%zsigi(jclass)
2283           zbnd2 = sec%zsigi(jclass+1)
2284        ENDIF
2285        !potential density classes transports
2286        IF( ( sec%zsigp(jclass)   .NE. 99._wp ) .AND. &
2287            ( sec%zsigp(jclass+1) .NE. 99._wp )       )THEN
2288           classe = 'DP      '
2289           zbnd1 = sec%zsigp(jclass)
2290           zbnd2 = sec%zsigp(jclass+1)
2291        ENDIF
2292        !depth classes transports
2293        IF( ( sec%zlay(jclass)    .NE. 99._wp ) .AND. &
2294            ( sec%zlay(jclass+1)  .NE. 99._wp )       )THEN
2295           classe = 'Z       '
2296           zbnd1 = sec%zlay(jclass)
2297           zbnd2 = sec%zlay(jclass+1)
2298        ENDIF
2299        !salinity classes transports
2300        IF( ( sec%zsal(jclass) .NE. 99._wp    ) .AND. &
2301            ( sec%zsal(jclass+1) .NE. 99._wp  )       )THEN
2302           classe = 'S       '
2303           zbnd1 = sec%zsal(jclass)
2304           zbnd2 = sec%zsal(jclass+1)   
2305        ENDIF
2306        !temperature classes transports
2307        IF( ( sec%ztem(jclass) .NE. 99._wp     ) .AND. &
2308            ( sec%ztem(jclass+1) .NE. 99._wp     )       ) THEN
2309           classe = 'T       '
2310           zbnd1 = sec%ztem(jclass)
2311           zbnd2 = sec%ztem(jclass+1)
2312        ENDIF
2313                 
2314        !write volume transport per class
2315        IF( lwp ) THEN
2316           
2317            IF  ( ln_dct_ascii ) THEN
2318                CALL FLUSH(numdct_NOOS)
2319
2320                !WRITE(numdct_NOOS,'(9e12.4E2)')   sec%transport( 1,jclass)+sec%transport( 2,jclass) , sec%transport( 1,jclass), sec%transport( 2,jclass), &
2321                !                                  sec%transport( 7,jclass)+sec%transport( 8,jclass) , sec%transport( 7,jclass), sec%transport( 8,jclass), &
2322                !                                  sec%transport( 9,jclass)+sec%transport(10,jclass) , sec%transport( 9,jclass), sec%transport(10,jclass)
2323                WRITE(numdct_NOOS,'(3F18.3,6e16.8E2)')   sec%transport( 1,jclass)+sec%transport( 2,jclass) , sec%transport( 1,jclass), sec%transport( 2,jclass), &
2324                                                         sec%transport( 7,jclass)+sec%transport( 8,jclass) , sec%transport( 7,jclass), sec%transport( 8,jclass), &
2325                                                         sec%transport( 9,jclass)+sec%transport(10,jclass) , sec%transport( 9,jclass), sec%transport(10,jclass)
2326            ELSE
2327
2328                CALL FLUSH(numdct_NOOS)
2329                WRITE(numdct_NOOS)   sec%transport( 1,jclass)+sec%transport( 2,jclass) , sec%transport( 1,jclass), sec%transport( 2,jclass), &
2330                                     sec%transport( 7,jclass)+sec%transport( 8,jclass) , sec%transport( 7,jclass), sec%transport( 8,jclass), &
2331                                     sec%transport( 9,jclass)+sec%transport(10,jclass) , sec%transport( 9,jclass), sec%transport(10,jclass)
2332            ENDIF
2333        ENDIF
2334
2335     ENDDO
2336     
2337     !IF  ( ln_dct_ascii ) THEN
2338        if ( lwp ) CALL FLUSH(numdct_NOOS)
2339     !ENDIF
2340
2341     CALL wrk_dealloc(nb_type_class , zsumclasses ) 
2342
2343  END SUBROUTINE dia_dct_wri_NOOS
2344
2345  SUBROUTINE dia_dct_wri_NOOS_h(hr,ksec,sec)
2346     !!-------------------------------------------------------------
2347     !! As routine dia_dct_wri_NOOS but for hourly output files
2348     !!
2349     !! Write transport output in numdct using NOOS formatting
2350     !!
2351     !! Purpose: Write  transports in ascii files
2352     !!
2353     !! Method:
2354     !!        1. Write volume transports in "volume_transport"
2355     !!           Unit: Sv : area * Velocity / 1.e6
2356     !!
2357     !!-------------------------------------------------------------
2358     !!arguments
2359     INTEGER, INTENT(IN)          :: hr          ! hour => effectively kt/12
2360     TYPE(SECTION), INTENT(INOUT) :: sec         ! section to write   
2361     INTEGER ,INTENT(IN)          :: ksec        ! section number
2362
2363     !!local declarations
2364     INTEGER               :: jclass,jhr            ! Dummy loop
2365     CHARACTER(len=2)      :: classe             ! Classname
2366     REAL(wp)              :: zbnd1,zbnd2        ! Class bounds
2367     REAL(wp)              :: zslope             ! section's slope coeff
2368     !
2369     REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace
2370     CHARACTER(len=3)      :: noos_sect_name             ! Classname
2371     CHARACTER(len=25)      :: noos_var_sect_name
2372     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   noos_iom_dummy
2373     INTEGER               :: IERR
2374     LOGICAL       ::   verbose     
2375     verbose = .false.
2376     
2377     !!-------------------------------------------------------------
2378
2379     IF( lwp .AND. verbose ) THEN
2380        WRITE(numout,*) " "
2381        WRITE(numout,*) "dia_dct_wri_NOOS_h: write transports through section Transect:",ksec-1," at timestep: ", hr
2382        WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
2383     ENDIF   
2384     
2385     CALL wrk_alloc(nb_type_class , zsumclasses ) 
2386     
2387     
2388     write (noos_sect_name, "(I03)")  ksec
2389     
2390     ALLOCATE( noos_iom_dummy(jpi,jpj,3),  STAT= ierr )
2391        IF( ierr /= 0 )   CALL ctl_stop( 'dia_dct_wri_NOOS_h: failed to allocate noos_iom_dummy array' )
2392
2393
2394
2395     
2396
2397     zsumclasses(:)=0._wp
2398     zslope = sec%slopeSection       
2399
2400     ! Sum up all classes, to give the total per type (pos vol trans, neg vol trans etc...)
2401     DO jclass=1,MAX(1,sec%nb_class-1)
2402        zsumclasses(1:nb_type_class)=zsumclasses(1:nb_type_class)+sec%transport_h(1:nb_type_class,jclass)
2403     ENDDO
2404       
2405     
2406     ! JT I think changing the sign of output according to the zslope is redundant
2407     
2408     !write volume transport per class
2409     ! Sum positive and vol trans for all classes in first cell of array
2410
2411     z_hr_output(ksec,1,1)= (zsumclasses(1)+zsumclasses(2))
2412     z_hr_output(ksec,2,1)= zsumclasses(1)
2413     z_hr_output(ksec,3,1)= zsumclasses(2)
2414
2415     ! Sum positive and vol trans for each classes in following cell of array
2416     DO jclass=1,MAX(1,sec%nb_class-1)
2417        z_hr_output(ksec,1,jclass+1)= (sec%transport_h(1,jclass)+sec%transport_h(2,jclass))
2418        z_hr_output(ksec,2,jclass+1)= sec%transport_h(1,jclass)
2419        z_hr_output(ksec,3,jclass+1)= sec%transport_h(2,jclass)
2420     ENDDO
2421
2422   
2423    IF( lwp )  THEN
2424        ! JT IF ( hr .eq. 48._wp ) THEN
2425        ! 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
2426        ! JT    DO jhr=25,48
2427        ! 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) )
2428        ! JT    ENDDO
2429        ! JT ENDIF
2430
2431
2432
2433        IF ( ln_dct_ascii ) THEN
2434            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
2435            WRITE(numdct_NOOS_h,'(11F18.3)')  z_hr_output(ksec,1,1), (z_hr_output(ksec,1,jclass+1), jclass=1,MAX(1,10) )
2436            WRITE(numdct_NOOS_h,'(11F18.3)')  z_hr_output(ksec,2,1), (z_hr_output(ksec,2,jclass+1), jclass=1,MAX(1,10) )
2437            WRITE(numdct_NOOS_h,'(11F18.3)')  z_hr_output(ksec,3,1), (z_hr_output(ksec,3,jclass+1), jclass=1,MAX(1,10) )
2438            CALL FLUSH(numdct_NOOS_h)
2439        ELSE
2440            WRITE(numdct_NOOS_h) nyear,nmonth,nday,MOD(hr,24),ksec-1,sec%nb_class-1
2441            WRITE(numdct_NOOS_h)  z_hr_output(ksec,1,1), (z_hr_output(ksec,1,jclass+1), jclass=1,MAX(1,10) )
2442            WRITE(numdct_NOOS_h)  z_hr_output(ksec,2,1), (z_hr_output(ksec,2,jclass+1), jclass=1,MAX(1,10) )
2443            WRITE(numdct_NOOS_h)  z_hr_output(ksec,3,1), (z_hr_output(ksec,3,jclass+1), jclass=1,MAX(1,10) )
2444            CALL FLUSH(numdct_NOOS_h)
2445        ENDIF
2446
2447
2448     ENDIF
2449
2450
2451     CALL wrk_dealloc(nb_type_class , zsumclasses )
2452     
2453     DEALLOCATE(noos_iom_dummy)
2454
2455
2456
2457  END SUBROUTINE dia_dct_wri_NOOS_h
2458
2459  SUBROUTINE dia_dct_wri(kt,ksec,sec)
2460     !!-------------------------------------------------------------
2461     !! Write transport output in numdct
2462     !!
2463     !! Purpose: Write  transports in ascii files
2464     !!
2465     !! Method:
2466     !!        1. Write volume transports in "volume_transport"
2467     !!           Unit: Sv : area * Velocity / 1.e6
2468     !!
2469     !!        2. Write heat transports in "heat_transport"
2470     !!           Unit: Peta W : area * Velocity * T * rhau * Cp / 1.e15
2471     !!
2472     !!        3. Write salt transports in "salt_transport"
2473     !!           Unit: 10^9 g m^3 / s : area * Velocity * S / 1.e6
2474     !!
2475     !!-------------------------------------------------------------
2476     !!arguments
2477     INTEGER, INTENT(IN)          :: kt          ! time-step
2478     TYPE(SECTION), INTENT(INOUT) :: sec         ! section to write   
2479     INTEGER ,INTENT(IN)          :: ksec        ! section number
2480
2481     !!local declarations
2482     INTEGER               :: jclass             ! Dummy loop
2483     CHARACTER(len=2)      :: classe             ! Classname
2484     REAL(wp)              :: zbnd1,zbnd2        ! Class bounds
2485     REAL(wp)              :: zslope             ! section's slope coeff
2486     !
2487     REAL(wp), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace
2488     !!-------------------------------------------------------------
2489     CALL wrk_alloc(nb_type_class , zsumclasses ) 
2490
2491     zsumclasses(:)=0._wp
2492     zslope = sec%slopeSection       
2493
2494 
2495     DO jclass=1,MAX(1,sec%nb_class-1)
2496
2497        classe   = 'N       '
2498        zbnd1   = 0._wp
2499        zbnd2   = 0._wp
2500        zsumclasses(1:nb_type_class)=zsumclasses(1:nb_type_class)+sec%transport(1:nb_type_class,jclass)
2501
2502   
2503        !insitu density classes transports
2504        IF( ( sec%zsigi(jclass)   .NE. 99._wp ) .AND. &
2505            ( sec%zsigi(jclass+1) .NE. 99._wp )       )THEN
2506           classe = 'DI       '
2507           zbnd1 = sec%zsigi(jclass)
2508           zbnd2 = sec%zsigi(jclass+1)
2509        ENDIF
2510        !potential density classes transports
2511        IF( ( sec%zsigp(jclass)   .NE. 99._wp ) .AND. &
2512            ( sec%zsigp(jclass+1) .NE. 99._wp )       )THEN
2513           classe = 'DP      '
2514           zbnd1 = sec%zsigp(jclass)
2515           zbnd2 = sec%zsigp(jclass+1)
2516        ENDIF
2517        !depth classes transports
2518        IF( ( sec%zlay(jclass)    .NE. 99._wp ) .AND. &
2519            ( sec%zlay(jclass+1)  .NE. 99._wp )       )THEN
2520           classe = 'Z       '
2521           zbnd1 = sec%zlay(jclass)
2522           zbnd2 = sec%zlay(jclass+1)
2523        ENDIF
2524        !salinity classes transports
2525        IF( ( sec%zsal(jclass) .NE. 99._wp    ) .AND. &
2526            ( sec%zsal(jclass+1) .NE. 99._wp  )       )THEN
2527           classe = 'S       '
2528           zbnd1 = sec%zsal(jclass)
2529           zbnd2 = sec%zsal(jclass+1)   
2530        ENDIF
2531        !temperature classes transports
2532        IF( ( sec%ztem(jclass) .NE. 99._wp     ) .AND. &
2533            ( sec%ztem(jclass+1) .NE. 99._wp     )       ) THEN
2534           classe = 'T       '
2535           zbnd1 = sec%ztem(jclass)
2536           zbnd2 = sec%ztem(jclass+1)
2537        ENDIF
2538                 
2539        !write volume transport per class
2540        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, &
2541                              jclass,classe,zbnd1,zbnd2,&
2542                              sec%transport(1,jclass),sec%transport(2,jclass), &
2543                              sec%transport(1,jclass)+sec%transport(2,jclass)
2544
2545        IF( sec%llstrpond )THEN
2546
2547           !write heat transport per class:
2548           WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope,  &
2549                              jclass,classe,zbnd1,zbnd2,&
2550                              sec%transport(7,jclass)*1000._wp*rcp/1.e15,sec%transport(8,jclass)*1000._wp*rcp/1.e15, &
2551                              ( sec%transport(7,jclass)+sec%transport(8,jclass) )*1000._wp*rcp/1.e15
2552           !write salt transport per class
2553           WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope,  &
2554                              jclass,classe,zbnd1,zbnd2,&
2555                              sec%transport(9,jclass)*1000._wp/1.e9,sec%transport(10,jclass)*1000._wp/1.e9,&
2556                              (sec%transport(9,jclass)+sec%transport(10,jclass))*1000._wp/1.e9
2557        ENDIF
2558
2559     ENDDO
2560
2561     zbnd1 = 0._wp
2562     zbnd2 = 0._wp
2563     jclass=0
2564
2565     !write total volume transport
2566     WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, &
2567                           jclass,"total",zbnd1,zbnd2,&
2568                           zsumclasses(1),zsumclasses(2),zsumclasses(1)+zsumclasses(2)
2569
2570     IF( sec%llstrpond )THEN
2571
2572        !write total heat transport
2573        WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope, &
2574                           jclass,"total",zbnd1,zbnd2,&
2575                           zsumclasses(7)* 1000._wp*rcp/1.e15,zsumclasses(8)* 1000._wp*rcp/1.e15,&
2576                           (zsumclasses(7)+zsumclasses(8) )* 1000._wp*rcp/1.e15
2577        !write total salt transport
2578        WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope, &
2579                           jclass,"total",zbnd1,zbnd2,&
2580                           zsumclasses(9)*1000._wp/1.e9,zsumclasses(10)*1000._wp/1.e9,&
2581                           (zsumclasses(9)+zsumclasses(10))*1000._wp/1.e9
2582     ENDIF
2583
2584     
2585     IF ( sec%ll_ice_section) THEN
2586        !write total ice volume transport
2587        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,&
2588                              jclass,"ice_vol",zbnd1,zbnd2,&
2589                              sec%transport(11,1),sec%transport(12,1),&
2590                              sec%transport(11,1)+sec%transport(12,1)
2591        !write total ice surface transport
2592        WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,&
2593                              jclass,"ice_surf",zbnd1,zbnd2,&
2594                              sec%transport(13,1),sec%transport(14,1), &
2595                              sec%transport(13,1)+sec%transport(14,1) 
2596     ENDIF
2597                                             
2598118 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3F12.4)
2599119 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6)
2600
2601     CALL wrk_dealloc(nb_type_class , zsumclasses ) 
2602  END SUBROUTINE dia_dct_wri
2603
2604   PURE  FUNCTION interp(ki, kj, kk, cd_point,var) 
2605  !!----------------------------------------------------------------------
2606  !!
2607  !!   Purpose: compute temperature/salinity/density at U-point or V-point
2608  !!   --------
2609  !!
2610  !!   Method:
2611  !!   ------
2612  !!
2613  !!   ====> full step and partial step
2614  !!
2615  !!
2616  !!    |    I          |    I+1           |    Z=temperature/salinity/density at U-poinT
2617  !!    |               |                  |
2618  !!  ----------------------------------------  1. Veritcale interpolation: compute zbis
2619  !!    |               |                  |       interpolation between ptab(I,J,K) and ptab(I,J,K+1)
2620  !!    |               |                  |       zbis =
2621  !!    |               |                  |      [ e3w(I+1,J,K)*ptab(I,J,K) + ( e3w(I,J,K) - e3w(I+1,J,K) ) * ptab(I,J,K-1) ]
2622  !!    |               |                  |      /[ e3w(I+1,J,K) + e3w(I,J,K) - e3w(I+1,J,K) ]
2623  !!    |               |                  |
2624  !!    |               |                  |    2. Horizontal interpolation: compute value at U/V point
2625  !!K-1 | ptab(I,J,K-1) |                  |       interpolation between zbis and ptab(I+1,J,K) 
2626  !!    |     .         |                  |
2627  !!    |     .         |                  |       interp = ( 0.5*zet2*zbis + 0.5*zet1*ptab(I+1,J,K) )/(0.5*zet2+0.5*zet1)
2628  !!    |     .         |                  |
2629  !!  ------------------------------------------
2630  !!    |     .         |                  |
2631  !!    |     .         |                  |
2632  !!    |     .         |                  |
2633  !!K   |    zbis.......U...ptab(I+1,J,K)  |
2634  !!    |     .         |                  |
2635  !!    | ptab(I,J,K)   |                  |
2636  !!    |               |------------------|
2637  !!    |               | partials         |
2638  !!    |               |  steps           |
2639  !!  -------------------------------------------
2640  !!    <----zet1------><----zet2--------->
2641  !!
2642  !!
2643  !!   ====>  s-coordinate
2644  !!     
2645  !!    |                |                  |   1. Compute distance between T1 and U points: SQRT( zdep1^2 + (0.5 * zet1 )^2
2646  !!    |                |                  |      Compute distance between T2 and U points: SQRT( zdep2^2 + (0.5 * zet2 )^2
2647  !!    |                | ptab(I+1,J,K)    |
2648  !!    |                |      T2          |   2. Interpolation between  T1 and T2 values at U point
2649  !!    |                |      ^           |   
2650  !!    |                |      | zdep2     |   
2651  !!    |                |      |           |   
2652  !!    |       ^        U      v           |
2653  !!    |       |        |                  |
2654  !!    |       | zdep1  |                  |   
2655  !!    |       v        |                  |
2656  !!    |      T1        |                  |
2657  !!    | ptab(I,J,K)    |                  |
2658  !!    |                |                  |
2659  !!    |                |                  |
2660  !!
2661  !!    <----zet1--------><----zet2--------->
2662  !!
2663  !!----------------------------------------------------------------------
2664  !*arguments
2665  INTEGER, INTENT(IN)                          :: ki, kj, kk   ! coordinate of point
2666  INTEGER, INTENT(IN)                          :: var   !  which variable
2667  CHARACTER(len=1), INTENT(IN)                 :: cd_point     ! type of point (U, V)
2668  REAL(wp)                                     :: interp       ! interpolated variable
2669
2670  !*local declations
2671  INTEGER :: ii1, ij1, ii2, ij2                                ! local integer
2672  REAL(wp):: ze3t, zfse3, zwgt1, zwgt2, zbis, zdepu            ! local real
2673  REAL(wp):: zet1, zet2                                        ! weight for interpolation
2674  REAL(wp):: zdep1,zdep2                                       ! differences of depth
2675  REAL(wp):: zmsk                                              ! mask value
2676  !!----------------------------------------------------------------------
2677
2678 
2679
2680  IF( cd_point=='U' )THEN
2681     ii1 = ki    ; ij1 = kj 
2682     ii2 = ki+1  ; ij2 = kj 
2683
2684     zet1=e1t(ii1,ij1)
2685     zet2=e1t(ii2,ij2)
2686     zmsk=umask(ii1,ij1,kk)
2687 
2688
2689  ELSE ! cd_point=='V'
2690     ii1 = ki    ; ij1 = kj 
2691     ii2 = ki    ; ij2 = kj+1 
2692
2693     zet1=e2t(ii1,ij1)
2694     zet2=e2t(ii2,ij2)
2695     zmsk=vmask(ii1,ij1,kk)
2696
2697  ENDIF
2698
2699  IF( ln_sco )THEN   ! s-coordinate case
2700
2701     zdepu = ( fsdept(ii1,ij1,kk) +  fsdept(ii2,ij2,kk) ) /2 
2702     zdep1 = fsdept(ii1,ij1,kk) - zdepu
2703     zdep2 = fsdept(ii2,ij2,kk) - zdepu
2704
2705     ! weights
2706     zwgt1 = SQRT( ( 0.5 * zet1 ) * ( 0.5 * zet1 ) + ( zdep1 * zdep1 ) )
2707     zwgt2 = SQRT( ( 0.5 * zet2 ) * ( 0.5 * zet2 ) + ( zdep2 * zdep2 ) )
2708 
2709     ! result
2710           SELECT CASE( var )
2711              CASE(0)  ;    interp = zmsk * ( zwgt2 *  tsn(ii1,ij1,kk,jp_tem) + zwgt1 *  tsn(ii1,ij1,kk,jp_tem) ) / ( zwgt2 + zwgt1 )
2712              CASE(1)  ;    interp = zmsk * ( zwgt2 *  tsn(ii1,ij1,kk,jp_sal) + zwgt1 *  tsn(ii1,ij1,kk,jp_sal) ) / ( zwgt2 + zwgt1 )
2713              CASE(2)  ;    interp = zmsk * ( zwgt2 *  rhop(ii1,ij1,kk) + zwgt1 *  rhop(ii1,ij1,kk) ) / ( zwgt2 + zwgt1 )
2714              CASE(3)  ;    interp = zmsk * ( zwgt2 *  (rhd(ii1,ij1,kk)*rau0+rau0) + zwgt1 *  (rhd(ii1,ij1,kk)*rau0+rau0) ) / ( zwgt2 + zwgt1 )
2715           END SELECT
2716
2717  ELSE       ! full step or partial step case
2718
2719#if defined key_vvl
2720
2721     ze3t  = fse3t_n(ii2,ij2,kk) - fse3t_n(ii1,ij1,kk) 
2722     zwgt1 = ( fse3w_n(ii2,ij2,kk) - fse3w_n(ii1,ij1,kk) ) / fse3w_n(ii2,ij2,kk)
2723     zwgt2 = ( fse3w_n(ii1,ij1,kk) - fse3w_n(ii2,ij2,kk) ) / fse3w_n(ii1,ij1,kk)
2724
2725#else
2726
2727     ze3t  = fse3t(ii2,ij2,kk)   - fse3t(ii1,ij1,kk) 
2728     zwgt1 = ( fse3w(ii2,ij2,kk) - fse3w(ii1,ij1,kk) ) / fse3w(ii2,ij2,kk)
2729     zwgt2 = ( fse3w(ii1,ij1,kk) - fse3w(ii2,ij2,kk) ) / fse3w(ii1,ij1,kk)
2730
2731#endif
2732
2733     IF(kk .NE. 1)THEN
2734
2735        IF( ze3t >= 0. )THEN 
2736           ! zbis
2737           SELECT CASE( var )
2738           CASE(0) 
2739                     zbis   =  tsn(ii2,ij2,kk,jp_tem) + zwgt1 *  (tsn(ii2,ij2,kk-1,jp_tem)-tsn(ii2,ij2,kk,jp_tem)   )
2740                     interp =  zmsk * ( zet2 * tsn(ii1,ij1,kk,jp_tem) + zet1 * zbis )/( zet1 + zet2 )
2741           CASE(1) 
2742                     zbis   =  tsn(ii2,ij2,kk,jp_sal) + zwgt1 *  (tsn(ii2,ij2,kk-1,jp_sal)-tsn(ii2,ij2,kk,jp_sal)   )
2743                     interp =  zmsk * ( zet2 * tsn(ii1,ij1,kk,jp_sal) + zet1 * zbis )/( zet1 + zet2 )
2744           CASE(2) 
2745                     zbis   =  rhop(ii2,ij2,kk) + zwgt1 *  (rhop(ii2,ij2,kk-1)-rhop(ii2,ij2,kk)   )
2746                     interp =  zmsk * ( zet2 * rhop(ii1,ij1,kk) + zet1 * zbis )/( zet1 + zet2 )
2747           CASE(3) 
2748                     zbis   =  (rhd(ii2,ij2,kk)*rau0+rau0) + zwgt1 *  ( (rhd(ii2,ij2,kk-1)*rau0+rau0)-(rhd(ii2,ij2,kk)*rau0+rau0)   )
2749                     interp =  zmsk * ( zet2 * (rhd(ii1,ij1,kk)*rau0+rau0) + zet1 * zbis )/( zet1 + zet2 )
2750           END SELECT
2751           ! result
2752        ELSE
2753           ! zbis
2754           SELECT CASE( var )
2755           CASE(0) 
2756                 zbis   = tsn(ii1,ij1,kk,jp_tem) + zwgt2 * ( tsn(ii1,ij1,kk-1,jp_tem) - tsn(ii1,ij2,kk,jp_tem) )
2757                 interp = zmsk * ( zet2 * zbis + zet1 * tsn(ii2,ij2,kk,jp_tem) )/( zet1 + zet2 )
2758           CASE(1) 
2759                 zbis   = tsn(ii1,ij1,kk,jp_sal) + zwgt2 * ( tsn(ii1,ij1,kk-1,jp_sal) - tsn(ii1,ij2,kk,jp_sal) )
2760                 interp = zmsk * ( zet2 * zbis + zet1 * tsn(ii2,ij2,kk,jp_sal) )/( zet1 + zet2 )
2761           CASE(2) 
2762                 zbis   = rhop(ii1,ij1,kk) + zwgt2 * ( rhop(ii1,ij1,kk-1) - rhop(ii1,ij2,kk) )
2763                 interp = zmsk * ( zet2 * zbis + zet1 * rhop(ii2,ij2,kk) )/( zet1 + zet2 )
2764           CASE(3) 
2765                 zbis   = (rhd(ii1,ij1,kk)*rau0+rau0) + zwgt2 * ( (rhd(ii1,ij1,kk-1)*rau0+rau0) - (rhd(ii1,ij2,kk)*rau0+rau0) )
2766                 interp = zmsk * ( zet2 * zbis + zet1 * (rhd(ii2,ij2,kk)*rau0+rau0) )/( zet1 + zet2 )
2767           END SELECT
2768        ENDIF   
2769
2770     ELSE
2771        SELECT CASE( var )
2772        CASE(0) 
2773           interp = zmsk * (  zet2 * tsn(ii1,ij1,kk,jp_tem) + zet1 * tsn(ii2,ij2,kk,jp_tem) )/( zet1 + zet2 )
2774        CASE(1) 
2775           interp = zmsk * (  zet2 * tsn(ii1,ij1,kk,jp_sal) + zet1 * tsn(ii2,ij2,kk,jp_sal) )/( zet1 + zet2 )
2776        CASE(2) 
2777           interp = zmsk * (  zet2 * rhop(ii1,ij1,kk) + zet1 * rhop(ii2,ij2,kk) )/( zet1 + zet2 )
2778        CASE(3) 
2779           interp = zmsk * (  zet2 * (rhd(ii1,ij1,kk)*rau0+rau0) + zet1 * (rhd(ii2,ij2,kk)*rau0+rau0) )/( zet1 + zet2 )
2780        END SELECT
2781     ENDIF
2782
2783  ENDIF
2784
2785  END FUNCTION interp
2786
2787#else
2788   !!----------------------------------------------------------------------
2789   !!   Default option :                                       Dummy module
2790   !!----------------------------------------------------------------------
2791   LOGICAL, PUBLIC, PARAMETER ::   lk_diadct = .FALSE.    !: diamht flag
2792   PUBLIC 
2793   !! $Id$
2794CONTAINS
2795
2796   SUBROUTINE dia_dct_init          ! Dummy routine
2797      WRITE(*,*) 'dia_dct_init: You should not have seen this print! error?', kt
2798   END SUBROUTINE dia_dct_init
2799
2800   SUBROUTINE dia_dct( kt )         ! Dummy routine
2801      INTEGER, INTENT( in ) :: kt   ! ocean time-step index
2802      WRITE(*,*) 'dia_dct: You should not have seen this print! error?', kt
2803   END SUBROUTINE dia_dct
2804#endif
2805
2806END MODULE diadct
Note: See TracBrowser for help on using the repository browser.