Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

diadct.F90

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

Download in other formats: