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

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

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

Context Navigation

source: branches/UKMO/AMM15_v3_6_STABLE_package_withNOOS/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 9472

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

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

Download in other formats: