New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

diadct.F90 in branches/UKMO/CO6_shelfclimate_fabm_noos/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: branches/UKMO/CO6_shelfclimate_fabm_noos/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 11264

Last change on this file since 11264 was 11264, checked in by hadjt, 5 years ago
Updated the maximum size of the NOOS cross-section size (from 150-275) in DIA/diadct.F90
Property svn:executable set to ``*
File size: 133.4 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: