source: branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis2/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 11038

Last change on this file since 11038 was 11038, checked in by rrenshaw, 17 months ago

code fixes

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