Context Navigation

diadct.F90 @ 7613

Last change on this file since 7613 was 7613, checked in by hadjt, 8 years ago
Added namelist entries to effectively switch off region mean, pea and NOOS code. Also use correct nemo file open subroutines, rahter than low level open statements with fixed file ids (37 and 73)
Property svn:executable set to ``*
File size: 102.1 KB

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

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

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

Context Navigation

source: branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 7613

Download in other formats: