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 NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIA – NEMO

Context Navigation

source: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/DIA/diadct.F90 @ 10876

Last change on this file since 10876 was 10009, checked in by gm, 6 years ago
#1911 (ENHANCE-04): RK3 branch - step II.1 time-level dimension on ssh
Property svn:keywords set to `Id`
File size: 58.1 KB

Line
1	MODULE diadct
2	!!======================================================================
3	!! * MODULE diadct *
4	!! Ocean diagnostics: Compute the transport trough a sec.
5	!!======================================================================
6	!! History : OPA ! 02/1999 (Y Drillet) original code
7	!! ! 10/2001 (Y Drillet, R Bourdalle Badie)
8	!! NEMO 1.0 ! 10/2005 (M Laborie) F90
9	!! 3.0 ! 04/2007 (G Garric) Ice sections
10	!! - ! 04/2007 (C Bricaud) test on sec%nb_point, initialisation of ztransp1,ztransp2,...
11	!! 3.4 ! 09/2011 (C Bricaud)
12	!!----------------------------------------------------------------------
13	#if defined key_diadct
14	!!----------------------------------------------------------------------
15	!! 'key_diadct' :
16	!!----------------------------------------------------------------------
17	!!----------------------------------------------------------------------
18	!! dia_dct : Compute the transport through a sec.
19	!! dia_dct_init : Read namelist.
20	!! readsec : Read sections description and pathway
21	!! removepoints : Remove points which are common to 2 procs
22	!! transport : Compute transport for each sections
23	!! dia_dct_wri : Write tranports results in ascii files
24	!! interp : Compute temperature/salinity/density at U-point or V-point
25	!!
26	!!----------------------------------------------------------------------
27	USE oce ! ocean dynamics and tracers
28	USE dom_oce ! ocean space and time domain
29	USE phycst ! physical constants
30	USE in_out_manager ! I/O manager
31	USE daymod ! calendar
32	USE dianam ! build name of file
33	USE lib_mpp ! distributed memory computing library
34	#if defined key_si3
35	USE ice
36	#endif
37	USE domvvl
38	USE timing ! preformance summary
39
40	IMPLICIT NONE
41	PRIVATE
42
43	PUBLIC dia_dct ! routine called by step.F90
44	PUBLIC dia_dct_init ! routine called by opa.F90
45	PUBLIC diadct_alloc ! routine called by nemo_init in nemogcm.F90
46	PRIVATE readsec
47	PRIVATE removepoints
48	PRIVATE transport
49	PRIVATE dia_dct_wri
50
51	LOGICAL, PUBLIC, PARAMETER :: lk_diadct = .TRUE. !: model-data diagnostics flag
52
53	INTEGER :: nn_dct ! Frequency of computation
54	INTEGER :: nn_dctwri ! Frequency of output
55	INTEGER :: nn_secdebug ! Number of the section to debug
56
57	INTEGER, PARAMETER :: nb_class_max = 10
58	INTEGER, PARAMETER :: nb_sec_max = 150
59	INTEGER, PARAMETER :: nb_point_max = 2000
60	INTEGER, PARAMETER :: nb_type_class = 10
61	INTEGER, PARAMETER :: nb_3d_vars = 3
62	INTEGER, PARAMETER :: nb_2d_vars = 2
63	INTEGER :: nb_sec
64
65	TYPE POINT_SECTION
66	INTEGER :: I,J
67	END TYPE POINT_SECTION
68
69	TYPE COORD_SECTION
70	REAL(wp) :: lon,lat
71	END TYPE COORD_SECTION
72
73	TYPE SECTION
74	CHARACTER(len=60) :: name ! name of the sec
75	LOGICAL :: llstrpond ! true if you want the computation of salt and
76	! heat transports
77	LOGICAL :: ll_ice_section ! ice surface and ice volume computation
78	LOGICAL :: ll_date_line ! = T if the section crosses the date-line
79	TYPE(COORD_SECTION), DIMENSION(2) :: coordSec ! longitude and latitude of the extremities of the sec
80	INTEGER :: nb_class ! number of boundaries for density classes
81	INTEGER, DIMENSION(nb_point_max) :: direction ! vector direction of the point in the section
82	CHARACTER(len=40),DIMENSION(nb_class_max) :: classname ! characteristics of the class
83	REAL(wp), DIMENSION(nb_class_max) :: zsigi ,&! in-situ density classes (99 if you don't want)
84	zsigp ,&! potential density classes (99 if you don't want)
85	zsal ,&! salinity classes (99 if you don't want)
86	ztem ,&! temperature classes(99 if you don't want)
87	zlay ! level classes (99 if you don't want)
88	REAL(wp), DIMENSION(nb_type_class,nb_class_max) :: transport ! transport output
89	REAL(wp) :: slopeSection ! slope of the section
90	INTEGER :: nb_point ! number of points in the section
91	TYPE(POINT_SECTION),DIMENSION(nb_point_max) :: listPoint ! list of points in the sections
92	END TYPE SECTION
93
94	TYPE(SECTION),DIMENSION(nb_sec_max) :: secs ! Array of sections
95
96	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: transports_3d
97	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: transports_2d
98
99	!!----------------------------------------------------------------------
100	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
101	!! $Id$
102	!! Software governed by the CeCILL licence (./LICENSE)
103	!!----------------------------------------------------------------------
104	CONTAINS
105
106	INTEGER FUNCTION diadct_alloc()
107	!!----------------------------------------------------------------------
108	!! * FUNCTION diadct_alloc *
109	!!----------------------------------------------------------------------
110	INTEGER :: ierr(2)
111	!!----------------------------------------------------------------------
112
113	ALLOCATE(transports_3d(nb_3d_vars,nb_sec_max,nb_point_max,jpk), STAT=ierr(1) )
114	ALLOCATE(transports_2d(nb_2d_vars,nb_sec_max,nb_point_max) , STAT=ierr(2) )
115
116	diadct_alloc = MAXVAL( ierr )
117	IF( diadct_alloc /= 0 ) CALL ctl_warn('diadct_alloc: failed to allocate arrays')
118
119	END FUNCTION diadct_alloc
120
121
122	SUBROUTINE dia_dct_init
123	!!---------------------------------------------------------------------
124	!! * ROUTINE diadct *
125	!!
126	!! ** Purpose: Read the namelist parameters
127	!! Open output files
128	!!
129	!!---------------------------------------------------------------------
130	INTEGER :: ios ! Local integer output status for namelist read
131	!!
132	NAMELIST/namdct/nn_dct,nn_dctwri,nn_secdebug
133	!!---------------------------------------------------------------------
134
135	REWIND( numnam_ref ) ! Namelist namdct in reference namelist : Diagnostic: transport through sections
136	READ ( numnam_ref, namdct, IOSTAT = ios, ERR = 901)
137	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdct in reference namelist', lwp )
138
139	REWIND( numnam_cfg ) ! Namelist namdct in configuration namelist : Diagnostic: transport through sections
140	READ ( numnam_cfg, namdct, IOSTAT = ios, ERR = 902 )
141	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdct in configuration namelist', lwp )
142	IF(lwm) WRITE ( numond, namdct )
143
144	IF( lwp ) THEN
145	WRITE(numout,*) " "
146	WRITE(numout,*) "diadct_init: compute transports through sections "
147	WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
148	WRITE(numout,*) " Frequency of computation: nn_dct = ",nn_dct
149	WRITE(numout,*) " Frequency of write: nn_dctwri = ",nn_dctwri
150
151	IF ( nn_secdebug .GE. 1 .AND. nn_secdebug .LE. nb_sec_max )THEN
152	WRITE(numout,*)" Debug section number: ", nn_secdebug
153	ELSE IF ( nn_secdebug == 0 )THEN ; WRITE(numout,*)" No section to debug"
154	ELSE IF ( nn_secdebug == -1 )THEN ; WRITE(numout,*)" Debug all sections"
155	ELSE ; WRITE(numout,*)" Wrong value for nn_secdebug : ",nn_secdebug
156	ENDIF
157
158	IF(nn_dct .GE. nn_dctwri .AND. MOD(nn_dct,nn_dctwri) .NE. 0) &
159	& CALL ctl_stop( 'diadct: nn_dct should be smaller and a multiple of nn_dctwri' )
160
161	ENDIF
162
163	!Read section_ijglobal.diadct
164	CALL readsec
165
166	!open output file
167	IF( lwm ) THEN
168	CALL ctl_opn( numdct_vol, 'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
169	CALL ctl_opn( numdct_heat, 'heat_transport' , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
170	CALL ctl_opn( numdct_salt, 'salt_transport' , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
171	ENDIF
172
173	! Initialise arrays to zero
174	transports_3d(:,:,:,:)=0.0
175	transports_2d(:,:,:) =0.0
176	!
177	END SUBROUTINE dia_dct_init
178
179
180	SUBROUTINE dia_dct( kt )
181	!!---------------------------------------------------------------------
182	!! * ROUTINE diadct *
183	!!
184	!! Purpose :: Compute section transports and write it in numdct files
185	!!
186	!! Method :: All arrays initialised to zero in dct_init
187	!! Each nn_dct time step call subroutine 'transports' for
188	!! each section to sum the transports over each grid cell.
189	!! Each nn_dctwri time step:
190	!! Divide the arrays by the number of summations to gain
191	!! an average value
192	!! Call dia_dct_sum to sum relevant grid boxes to obtain
193	!! totals for each class (density, depth, temp or sal)
194	!! Call dia_dct_wri to write the transports into file
195	!! Reinitialise all relevant arrays to zero
196	!!---------------------------------------------------------------------
197	INTEGER, INTENT(in) :: kt
198	!
199	INTEGER :: jsec ! loop on sections
200	INTEGER :: itotal ! nb_sec_maxnb_type_classnb_class_max
201	LOGICAL :: lldebug =.FALSE. ! debug a section
202	INTEGER , DIMENSION(1) :: ish ! work array for mpp_sum
203	INTEGER , DIMENSION(3) :: ish2 ! "
204	REAL(wp), ALLOCATABLE, DIMENSION(:) :: zwork ! "
205	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:):: zsum ! "
206	!!---------------------------------------------------------------------
207	!
208	IF( ln_timing ) CALL timing_start('dia_dct')
209
210	IF( lk_mpp )THEN
211	itotal = nb_sec_maxnb_type_classnb_class_max
212	ALLOCATE( zwork(itotal) , zsum(nb_sec_max,nb_type_class,nb_class_max) )
213	ENDIF
214
215	! Initialise arrays
216	zwork(:) = 0.0
217	zsum(:,:,:) = 0.0
218
219	IF( lwp .AND. kt==nit000+nn_dct-1 ) THEN
220	WRITE(numout,*) " "
221	WRITE(numout,*) "diadct: compute transport"
222	WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~~~~~"
223	WRITE(numout,*) "nb_sec = ",nb_sec
224	ENDIF
225
226
227	! Compute transport and write only at nn_dctwri
228	IF( MOD(kt,nn_dct)==0 ) THEN
229
230	DO jsec=1,nb_sec
231
232	!debug this section computing ?
233	lldebug=.FALSE.
234	IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND. kt==nit000+nn_dct-1 ) lldebug=.TRUE.
235
236	!Compute transport through section
237	CALL transport(secs(jsec),lldebug,jsec)
238
239	ENDDO
240
241	IF( MOD(kt,nn_dctwri)==0 )THEN
242
243	IF( kt==nit000+nn_dctwri-1 )WRITE(numout,*)" diadct: average transports and write at kt = ",kt
244
245	!! divide arrays by nn_dctwri/nn_dct to obtain average
246	transports_3d(:,:,:,:)=transports_3d(:,:,:,:)/(nn_dctwri/nn_dct)
247	transports_2d(:,:,:) =transports_2d(:,:,:) /(nn_dctwri/nn_dct)
248
249	! Sum over each class
250	DO jsec=1,nb_sec
251	CALL dia_dct_sum(secs(jsec),jsec)
252	ENDDO
253
254	!Sum on all procs
255	IF( lk_mpp )THEN
256	ish(1) = nb_sec_maxnb_type_classnb_class_max
257	ish2 = (/nb_sec_max,nb_type_class,nb_class_max/)
258	DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport(:,:) ; ENDDO
259	zwork(:)= RESHAPE(zsum(:,:,:), ish )
260	CALL mpp_sum(zwork, ish(1))
261	zsum(:,:,:)= RESHAPE(zwork,ish2)
262	DO jsec=1,nb_sec ; secs(jsec)%transport(:,:) = zsum(jsec,:,:) ; ENDDO
263	ENDIF
264
265	!Write the transport
266	DO jsec=1,nb_sec
267
268	IF( lwm )CALL dia_dct_wri(kt,jsec,secs(jsec))
269
270	!nullify transports values after writing
271	transports_3d(:,jsec,:,:)=0.
272	transports_2d(:,jsec,: )=0.
273	secs(jsec)%transport(:,:)=0.
274
275	ENDDO
276
277	ENDIF
278
279	ENDIF
280
281	IF( lk_mpp )THEN
282	itotal = nb_sec_maxnb_type_classnb_class_max
283	DEALLOCATE( zwork , zsum )
284	ENDIF
285
286	IF( ln_timing ) CALL timing_stop('dia_dct')
287	!
288	END SUBROUTINE dia_dct
289
290
291	SUBROUTINE readsec
292	!!---------------------------------------------------------------------
293	!! * ROUTINE readsec *
294	!!
295	!! ** Purpose:
296	!! Read a binary file(section_ijglobal.diadct)
297	!! generated by the tools "NEMOGCM/TOOLS/SECTIONS_DIADCT"
298	!!
299	!!
300	!!---------------------------------------------------------------------
301	INTEGER :: iptglo , iptloc ! Global and local number of points for a section
302	INTEGER :: isec, iiglo, ijglo, iiloc, ijloc,iost,i1 ,i2 ! temporary integer
303	INTEGER :: jsec, jpt ! dummy loop indices
304	INTEGER, DIMENSION(2) :: icoord
305	LOGICAL :: llbon, lldebug ! local logical
306	CHARACTER(len=160) :: clname ! filename
307	CHARACTER(len=200) :: cltmp
308	CHARACTER(len=200) :: clformat !automatic format
309	TYPE(POINT_SECTION),DIMENSION(nb_point_max) ::coordtemp !contains listpoints coordinates read in the file
310	INTEGER, DIMENSION(nb_point_max) :: directemp !contains listpoints directions read in the files
311	!!-------------------------------------------------------------------------------------
312
313	!open input file
314	!---------------
315	CALL ctl_opn( numdct_in, 'section_ijglobal.diadct', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
316
317	!---------------
318	!Read input file
319	!---------------
320
321	DO jsec=1,nb_sec_max !loop on the nb_sec sections
322
323	IF ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) &
324	& WRITE(numout,*)'debuging for section number: ',jsec
325
326	!initialization
327	!---------------
328	secs(jsec)%name=''
329	secs(jsec)%llstrpond = .FALSE. ; secs(jsec)%ll_ice_section = .FALSE.
330	secs(jsec)%ll_date_line = .FALSE. ; secs(jsec)%nb_class = 0
331	secs(jsec)%zsigi = 99._wp ; secs(jsec)%zsigp = 99._wp
332	secs(jsec)%zsal = 99._wp ; secs(jsec)%ztem = 99._wp
333	secs(jsec)%zlay = 99._wp
334	secs(jsec)%transport = 0._wp ; secs(jsec)%nb_point = 0
335
336	!read section's number / name / computing choices / classes / slopeSection / points number
337	!-----------------------------------------------------------------------------------------
338	READ(numdct_in,iostat=iost)isec
339	IF (iost .NE. 0 )EXIT !end of file
340	WRITE(cltmp,'(a,i4.4,a,i4.4)')'diadct: read sections : Problem of section number: isec= ',isec,' and jsec= ',jsec
341	IF( jsec .NE. isec ) CALL ctl_stop( cltmp )
342
343	IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )WRITE(numout,*)"isec ",isec
344
345	READ(numdct_in)secs(jsec)%name
346	READ(numdct_in)secs(jsec)%llstrpond
347	READ(numdct_in)secs(jsec)%ll_ice_section
348	READ(numdct_in)secs(jsec)%ll_date_line
349	READ(numdct_in)secs(jsec)%coordSec
350	READ(numdct_in)secs(jsec)%nb_class
351	READ(numdct_in)secs(jsec)%zsigi
352	READ(numdct_in)secs(jsec)%zsigp
353	READ(numdct_in)secs(jsec)%zsal
354	READ(numdct_in)secs(jsec)%ztem
355	READ(numdct_in)secs(jsec)%zlay
356	READ(numdct_in)secs(jsec)%slopeSection
357	READ(numdct_in)iptglo
358
359	!debug
360	!-----
361
362	IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )THEN
363
364	WRITE(clformat,'(a,i2,a)') '(A40,', nb_class_max,'(f8.3,1X))'
365
366	WRITE(numout,*) " Section name : ",TRIM(secs(jsec)%name)
367	WRITE(numout,*) " Compute heat and salt transport ? ",secs(jsec)%llstrpond
368	WRITE(numout,*) " Compute ice transport ? ",secs(jsec)%ll_ice_section
369	WRITE(numout,*) " Section crosses date-line ? ",secs(jsec)%ll_date_line
370	WRITE(numout,*) " Slope section : ",secs(jsec)%slopeSection
371	WRITE(numout,*) " Number of points in the section: ",iptglo
372	WRITE(numout,*) " Number of classes ",secs(jsec)%nb_class
373	WRITE(numout,clformat)" Insitu density classes : ",secs(jsec)%zsigi
374	WRITE(numout,clformat)" Potential density classes : ",secs(jsec)%zsigp
375	WRITE(numout,clformat)" Salinity classes : ",secs(jsec)%zsal
376	WRITE(numout,clformat)" Temperature classes : ",secs(jsec)%ztem
377	WRITE(numout,clformat)" Depth classes : ",secs(jsec)%zlay
378	ENDIF
379
380	IF( iptglo /= 0 )THEN
381
382	!read points'coordinates and directions
383	!--------------------------------------
384	coordtemp(:) = POINT_SECTION(0,0) !list of points read
385	directemp(:) = 0 !value of directions of each points
386	DO jpt=1,iptglo
387	READ(numdct_in) i1, i2
388	coordtemp(jpt)%I = i1
389	coordtemp(jpt)%J = i2
390	ENDDO
391	READ(numdct_in) directemp(1:iptglo)
392
393	!debug
394	!-----
395	IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )THEN
396	WRITE(numout,*)" List of points in global domain:"
397	DO jpt=1,iptglo
398	WRITE(numout,*)' # I J ',jpt,coordtemp(jpt),directemp(jpt)
399	ENDDO
400	ENDIF
401
402	!Now each proc selects only points that are in its domain:
403	!--------------------------------------------------------
404	iptloc = 0 ! initialize number of points selected
405	DO jpt = 1, iptglo ! loop on listpoint read in the file
406	!
407	iiglo=coordtemp(jpt)%I ! global coordinates of the point
408	ijglo=coordtemp(jpt)%J ! "
409
410	IF( iiglo==jpiglo .AND. nimpp==1 ) iiglo = 2 !!gm BUG: Hard coded periodicity !
411
412	iiloc=iiglo-nimpp+1 ! local coordinates of the point
413	ijloc=ijglo-njmpp+1 ! "
414
415	!verify if the point is on the local domain:(1,nlei)*(1,nlej)
416	IF( iiloc >= 1 .AND. iiloc <= nlei .AND. &
417	ijloc >= 1 .AND. ijloc <= nlej )THEN
418	iptloc = iptloc + 1 ! count local points
419	secs(jsec)%listPoint(iptloc) = POINT_SECTION(mi0(iiglo),mj0(ijglo)) ! store local coordinates
420	secs(jsec)%direction(iptloc) = directemp(jpt) ! store local direction
421	ENDIF
422	!
423	END DO
424
425	secs(jsec)%nb_point=iptloc !store number of section's points
426
427	!debug
428	!-----
429	IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )THEN
430	WRITE(numout,*)" List of points selected by the proc:"
431	DO jpt = 1,iptloc
432	iiglo = secs(jsec)%listPoint(jpt)%I + nimpp - 1
433	ijglo = secs(jsec)%listPoint(jpt)%J + njmpp - 1
434	WRITE(numout,*)' # I J : ',iiglo,ijglo
435	ENDDO
436	ENDIF
437
438	IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN
439	DO jpt = 1,iptloc
440	iiglo = secs(jsec)%listPoint(jpt)%I + nimpp - 1
441	ijglo = secs(jsec)%listPoint(jpt)%J + njmpp - 1
442	ENDDO
443	ENDIF
444
445	!remove redundant points between processors
446	!------------------------------------------
447	lldebug = .FALSE. ; IF ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) lldebug = .TRUE.
448	IF( iptloc .NE. 0 )THEN
449	CALL removepoints(secs(jsec),'I','top_list',lldebug)
450	CALL removepoints(secs(jsec),'I','bot_list',lldebug)
451	CALL removepoints(secs(jsec),'J','top_list',lldebug)
452	CALL removepoints(secs(jsec),'J','bot_list',lldebug)
453	ENDIF
454	IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN
455	DO jpt = 1,secs(jsec)%nb_point
456	iiglo = secs(jsec)%listPoint(jpt)%I + nimpp - 1
457	ijglo = secs(jsec)%listPoint(jpt)%J + njmpp - 1
458	ENDDO
459	ENDIF
460
461	!debug
462	!-----
463	IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )THEN
464	WRITE(numout,*)" List of points after removepoints:"
465	iptloc = secs(jsec)%nb_point
466	DO jpt = 1,iptloc
467	iiglo = secs(jsec)%listPoint(jpt)%I + nimpp - 1
468	ijglo = secs(jsec)%listPoint(jpt)%J + njmpp - 1
469	WRITE(numout,*)' # I J : ',iiglo,ijglo
470	CALL FLUSH(numout)
471	ENDDO
472	ENDIF
473
474	ELSE ! iptglo = 0
475	IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )&
476	WRITE(numout,*)' No points for this section.'
477	ENDIF
478
479	ENDDO !end of the loop on jsec
480
481	nb_sec = jsec-1 !number of section read in the file
482	!
483	END SUBROUTINE readsec
484
485
486	SUBROUTINE removepoints(sec,cdind,cdextr,ld_debug)
487	!!---------------------------------------------------------------------------
488	!! *** function removepoints
489	!!
490	!! ** Purpose :: Remove points which are common to 2 procs
491	!!
492	!----------------------------------------------------------------------------
493	TYPE(SECTION),INTENT(INOUT) :: sec
494	CHARACTER(len=1),INTENT(IN) :: cdind ! = 'I'/'J'
495	CHARACTER(len=8),INTENT(IN) :: cdextr ! = 'top_list'/'bot_list'
496	LOGICAL,INTENT(IN) :: ld_debug
497	!
498	INTEGER :: iextr ,& !extremity of listpoint that we verify
499	iind ,& !coord of listpoint that we verify
500	itest ,& !indice value of the side of the domain
501	!where points could be redundant
502	isgn ,& ! isgn= 1 : scan listpoint from start to end
503	! isgn=-1 : scan listpoint from end to start
504	istart,iend !first and last points selected in listpoint
505	INTEGER :: jpoint !loop on list points
506	INTEGER, DIMENSION(nb_point_max) :: idirec !contains temporary sec%direction
507	INTEGER, DIMENSION(2,nb_point_max) :: icoord !contains temporary sec%listpoint
508	!----------------------------------------------------------------------------
509	!
510	IF( ld_debug )WRITE(numout,*)' -------------------------'
511	IF( ld_debug )WRITE(numout,*)' removepoints in listpoint'
512
513	!iextr=extremity of list_point that we verify
514	IF ( cdextr=='bot_list' )THEN ; iextr=1 ; isgn=1
515	ELSE IF ( cdextr=='top_list' )THEN ; iextr=sec%nb_point ; isgn=-1
516	ELSE ; CALL ctl_stop("removepoints :Wrong value for cdextr")
517	ENDIF
518
519	!which coordinate shall we verify ?
520	IF ( cdind=='I' )THEN ; itest=nlei ; iind=1
521	ELSE IF ( cdind=='J' )THEN ; itest=nlej ; iind=2
522	ELSE ; CALL ctl_stop("removepoints :Wrong value for cdind")
523	ENDIF
524
525	IF( ld_debug )THEN
526	WRITE(numout,*)' case: coord/list extr/domain side'
527	WRITE(numout,*)' ', cdind,' ',cdextr,' ',itest
528	WRITE(numout,*)' Actual number of points: ',sec%nb_point
529	ENDIF
530
531	icoord(1,1:nb_point_max) = sec%listPoint%I
532	icoord(2,1:nb_point_max) = sec%listPoint%J
533	idirec = sec%direction
534	sec%listPoint = POINT_SECTION(0,0)
535	sec%direction = 0
536
537	jpoint=iextr+isgn
538	DO WHILE( jpoint .GE. 1 .AND. jpoint .LE. sec%nb_point )
539	IF( icoord( iind,jpoint-isgn ) == itest .AND. icoord( iind,jpoint ) == itest )THEN ; jpoint=jpoint+isgn
540	ELSE ; EXIT
541	ENDIF
542	ENDDO
543
544	IF( cdextr=='bot_list')THEN ; istart=jpoint-1 ; iend=sec%nb_point
545	ELSE ; istart=1 ; iend=jpoint+1
546	ENDIF
547
548	sec%listPoint(1:1+iend-istart)%I = icoord(1,istart:iend)
549	sec%listPoint(1:1+iend-istart)%J = icoord(2,istart:iend)
550	sec%direction(1:1+iend-istart) = idirec(istart:iend)
551	sec%nb_point = iend-istart+1
552
553	IF( ld_debug )THEN
554	WRITE(numout,*)' Number of points after removepoints :',sec%nb_point
555	WRITE(numout,*)' sec%direction after removepoints :',sec%direction(1:sec%nb_point)
556	ENDIF
557	!
558	END SUBROUTINE removepoints
559
560
561	SUBROUTINE transport(sec,ld_debug,jsec)
562	!!-------------------------------------------------------------------------------------------
563	!! * ROUTINE transport *
564	!!
565	!! Purpose :: Compute the transport for each point in a section
566	!!
567	!! Method :: Loop over each segment, and each vertical level and add the transport
568	!! Be aware :
569	!! One section is a sum of segments
570	!! One segment is defined by 2 consecutive points in sec%listPoint
571	!! All points of sec%listPoint are positioned on the F-point of the cell
572	!!
573	!! There are two loops:
574	!! loop on the segment between 2 nodes
575	!! loop on the level jk !!
576	!!
577	!! Output :: Arrays containing the volume,density,heat,salt transports for each i
578	!! point in a section, summed over each nn_dct.
579	!!
580	!!-------------------------------------------------------------------------------------------
581	TYPE(SECTION),INTENT(INOUT) :: sec
582	LOGICAL ,INTENT(IN) :: ld_debug
583	INTEGER ,INTENT(IN) :: jsec ! numeric identifier of section
584	!
585	INTEGER :: jk, jseg, jclass,jl, isgnu, isgnv ! loop on level/segment/classes/ice categories
586	REAL(wp):: zumid, zvmid, zumid_ice, zvmid_ice ! U/V ocean & ice velocity on a cell segment
587	REAL(wp):: zTnorm ! transport of velocity through one cell's sides
588	REAL(wp):: ztn, zsn, zrhoi, zrhop, zssh, zdep ! temperature/salinity/potential density/ssh/depth at u/v point
589	TYPE(POINT_SECTION) :: k
590	!!--------------------------------------------------------
591	!
592	IF( ld_debug )WRITE(numout,*)' Compute transport'
593
594	!---------------------------!
595	! COMPUTE TRANSPORT !
596	!---------------------------!
597	IF(sec%nb_point .NE. 0)THEN
598
599	!----------------------------------------------------------------------------------------------------
600	!Compute sign for velocities:
601	!
602	!convention:
603	! non horizontal section: direction + is toward left hand of section
604	! horizontal section: direction + is toward north of section
605	!
606	!
607	! slopeSection < 0 slopeSection > 0 slopeSection=inf slopeSection=0
608	! ---------------- ----------------- --------------- --------------
609	!
610	! isgnv=1 direction +
611	! ______ _____ ______
612	! \| //\| \| \| direction +
613	! \| isgnu=1 // \| \|isgnu=1 \|isgnu=1 /\|\
614	! \|_______ // ______\| \\ \| ---\ \|
615	! \| \| isgnv=-1 \\ \| \| ---/ direction + ____________
616	! \| \| __\\\| \|
617	! \| \| direction + \| isgnv=1
618	!
619	!----------------------------------------------------------------------------------------------------
620	isgnu = 1
621	IF( sec%slopeSection .GT. 0 ) THEN ; isgnv = -1
622	ELSE ; isgnv = 1
623	ENDIF
624	IF( sec%slopeSection .GE. 9999. ) isgnv = 1
625
626	IF( ld_debug )write(numout,*)"sec%slopeSection isgnu isgnv ",sec%slopeSection,isgnu,isgnv
627
628	!--------------------------------------!
629	! LOOP ON THE SEGMENT BETWEEN 2 NODES !
630	!--------------------------------------!
631	DO jseg=1,MAX(sec%nb_point-1,0)
632
633	!-------------------------------------------------------------------------------------------
634	! Select the appropriate coordinate for computing the velocity of the segment
635	!
636	! CASE(0) Case (2)
637	! ------- --------
638	! listPoint(jseg) listPoint(jseg+1) listPoint(jseg) F(i,j)
639	! F(i,j)----------V(i+1,j)-------F(i+1,j) \|
640	! \|
641	! \|
642	! \|
643	! Case (3) U(i,j)
644	! -------- \|
645	! \|
646	! listPoint(jseg+1) F(i,j+1) \|
647	! \| \|
648	! \| \|
649	! \| listPoint(jseg+1) F(i,j-1)
650	! \|
651	! \|
652	! U(i,j+1)
653	! \| Case(1)
654	! \| ------
655	! \|
656	! \| listPoint(jseg+1) listPoint(jseg)
657	! \| F(i-1,j)-----------V(i,j) -------f(jseg)
658	! listPoint(jseg) F(i,j)
659	!
660	!-------------------------------------------------------------------------------------------
661
662	SELECT CASE( sec%direction(jseg) )
663	CASE(0) ; k = sec%listPoint(jseg)
664	CASE(1) ; k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
665	CASE(2) ; k = sec%listPoint(jseg)
666	CASE(3) ; k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
667	END SELECT
668
669	!---------------------------\|
670	! LOOP ON THE LEVEL \|
671	!---------------------------\|
672	DO jk = 1, mbkt(k%I,k%J) !Sum of the transport on the vertical
673	! ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
674	SELECT CASE( sec%direction(jseg) )
675	CASE(0,1)
676	ztn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) )
677	zsn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) )
678	zrhop = interp(k%I,k%J,jk,'V',rhop)
679	zrhoi = interp(k%I,k%J,jk,'V',rhd*rho0+rho0)
680	zssh = 0.5( ssh(k%I,k%J,Nnn) + ssh(k%I,k%J+1,Nnn) ) vmask(k%I,k%J,1)
681	CASE(2,3)
682	ztn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) )
683	zsn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) )
684	zrhop = interp(k%I,k%J,jk,'U',rhop)
685	zrhoi = interp(k%I,k%J,jk,'U',rhd*rho0+rho0)
686	zssh = 0.5( ssh(k%I,k%J,Nnn) + ssh(k%I+1,k%J,Nnn) ) umask(k%I,k%J,1)
687	END SELECT
688	!
689	zdep= gdept_n(k%I,k%J,jk)
690
691	SELECT CASE( sec%direction(jseg) ) !compute velocity with the correct direction
692	CASE(0,1)
693	zumid=0._wp
694	zvmid=isgnvvn(k%I,k%J,jk)vmask(k%I,k%J,jk)
695	CASE(2,3)
696	zumid=isgnuun(k%I,k%J,jk)umask(k%I,k%J,jk)
697	zvmid=0._wp
698	END SELECT
699
700	!zTnorm=transport through one cell;
701	!velocity* cell's length * cell's thickness
702	zTnorm = zumide2u(k%I,k%J) e3u_n(k%I,k%J,jk) &
703	& + zvmide1v(k%I,k%J) e3v_n(k%I,k%J,jk)
704
705	!!gm THIS is WRONG no transport due to ssh in linear free surface case !!!!!
706	! IF( ln_linssh ) THEN !add transport due to free surface
707	! IF( jk==1 ) THEN
708	! zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zssh * umask(k%I,k%J,jk) &
709	! & + zvmid* e1v(k%I,k%J) * zssh * vmask(k%I,k%J,jk)
710	! ENDIF
711	! ENDIF
712	!!gm end
713	!COMPUTE TRANSPORT
714
715	transports_3d(1,jsec,jseg,jk) = transports_3d(1,jsec,jseg,jk) + zTnorm
716
717	IF( sec%llstrpond ) THEN
718	transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk) + zTnorm * ztn * zrhop * rcp
719	transports_3d(3,jsec,jseg,jk) = transports_3d(3,jsec,jseg,jk) + zTnorm * zsn * zrhop * 0.001
720	ENDIF
721
722	END DO !end of loop on the level
723
724	#if defined key_si3
725
726	!ICE CASE
727	!------------
728	IF( sec%ll_ice_section )THEN
729	SELECT CASE (sec%direction(jseg))
730	CASE(0)
731	zumid_ice = 0
732	zvmid_ice = isgnv0.5(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
733	CASE(1)
734	zumid_ice = 0
735	zvmid_ice = isgnv0.5(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
736	CASE(2)
737	zvmid_ice = 0
738	zumid_ice = isgnu0.5(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
739	CASE(3)
740	zvmid_ice = 0
741	zumid_ice = isgnu0.5(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
742	END SELECT
743
744	zTnorm=zumid_icee2u(k%I,k%J)+zvmid_icee1v(k%I,k%J)
745
746	#if defined key_si3
747	DO jl=1,jpl
748	transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)* &
749	a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) * &
750	( h_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) + &
751	h_s(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) )
752
753	transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)* &
754	a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl)
755	END DO
756	#endif
757
758	ENDIF !end of ice case
759	#endif
760
761	END DO !end of loop on the segment
762
763	ENDIF !end of sec%nb_point =0 case
764	!
765	END SUBROUTINE transport
766
767
768	SUBROUTINE dia_dct_sum(sec,jsec)
769	!!-------------------------------------------------------------
770	!! Purpose: Average the transport over nn_dctwri time steps
771	!! and sum over the density/salinity/temperature/depth classes
772	!!
773	!! Method: Sum over relevant grid cells to obtain values
774	!! for each class
775	!! There are several loops:
776	!! loop on the segment between 2 nodes
777	!! loop on the level jk
778	!! loop on the density/temperature/salinity/level classes
779	!! test on the density/temperature/salinity/level
780	!!
781	!! Note: Transport through a given section is equal to the sum of transports
782	!! computed on each proc.
783	!! On each proc,transport is equal to the sum of transport computed through
784	!! segments linking each point of sec%listPoint with the next one.
785	!!
786	!!-------------------------------------------------------------
787	TYPE(SECTION),INTENT(INOUT) :: sec
788	INTEGER ,INTENT(IN) :: jsec ! numeric identifier of section
789
790	TYPE(POINT_SECTION) :: k
791	INTEGER :: jk,jseg,jclass ! dummy variables for looping on level/segment/classes
792	REAL(wp) :: ztn, zsn, zrhoi, zrhop, zssh, zdep ! temperature/salinity/ssh/potential density /depth at u/v point
793	!!-------------------------------------------------------------
794
795	!! Sum the relevant segments to obtain values for each class
796	IF(sec%nb_point .NE. 0)THEN
797
798	!--------------------------------------!
799	! LOOP ON THE SEGMENT BETWEEN 2 NODES !
800	!--------------------------------------!
801	DO jseg=1,MAX(sec%nb_point-1,0)
802
803	!-------------------------------------------------------------------------------------------
804	! Select the appropriate coordinate for computing the velocity of the segment
805	!
806	! CASE(0) Case (2)
807	! ------- --------
808	! listPoint(jseg) listPoint(jseg+1) listPoint(jseg) F(i,j)
809	! F(i,j)----------V(i+1,j)-------F(i+1,j) \|
810	! \|
811	! \|
812	! \|
813	! Case (3) U(i,j)
814	! -------- \|
815	! \|
816	! listPoint(jseg+1) F(i,j+1) \|
817	! \| \|
818	! \| \|
819	! \| listPoint(jseg+1) F(i,j-1)
820	! \|
821	! \|
822	! U(i,j+1)
823	! \| Case(1)
824	! \| ------
825	! \|
826	! \| listPoint(jseg+1) listPoint(jseg)
827	! \| F(i-1,j)-----------V(i,j) -------f(jseg)
828	! listPoint(jseg) F(i,j)
829	!
830	!-------------------------------------------------------------------------------------------
831
832	SELECT CASE( sec%direction(jseg) )
833	CASE(0) ; k = sec%listPoint(jseg)
834	CASE(1) ; k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
835	CASE(2) ; k = sec%listPoint(jseg)
836	CASE(3) ; k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
837	END SELECT
838
839	!---------------------------\|
840	! LOOP ON THE LEVEL \|
841	!---------------------------\|
842	!Sum of the transport on the vertical
843	DO jk=1,mbkt(k%I,k%J)
844
845	! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
846	SELECT CASE( sec%direction(jseg) )
847	CASE(0,1)
848	ztn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) )
849	zsn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) )
850	zrhop = interp(k%I,k%J,jk,'V',rhop)
851	zrhoi = interp(k%I,k%J,jk,'V',rhd*rho0+rho0)
852
853	CASE(2,3)
854	ztn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) )
855	zsn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) )
856	zrhop = interp(k%I,k%J,jk,'U',rhop)
857	zrhoi = interp(k%I,k%J,jk,'U',rhd*rho0+rho0)
858	zssh = 0.5( ssh(k%I,k%J,Nnn) + ssh(k%I+1,k%J,Nnn) ) umask(k%I,k%J,1)
859	END SELECT
860
861	zdep= gdept_n(k%I,k%J,jk)
862
863	!-------------------------------
864	! LOOP ON THE DENSITY CLASSES \|
865	!-------------------------------
866	!The computation is made for each density/temperature/salinity/depth class
867	DO jclass=1,MAX(1,sec%nb_class-1)
868
869	!----------------------------------------------!
870	!TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL!
871	!----------------------------------------------!
872
873	IF ( ( &
874	((( zrhop .GE. (sec%zsigp(jclass)+1000. )) .AND. &
875	( zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR. &
876	( sec%zsigp(jclass) .EQ. 99.)) .AND. &
877
878	((( zrhoi .GE. (sec%zsigi(jclass) + 1000. )) .AND. &
879	( zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR. &
880	( sec%zsigi(jclass) .EQ. 99.)) .AND. &
881
882	((( zsn .GT. sec%zsal(jclass)) .AND. &
883	( zsn .LE. sec%zsal(jclass+1))) .OR. &
884	( sec%zsal(jclass) .EQ. 99.)) .AND. &
885
886	((( ztn .GE. sec%ztem(jclass)) .AND. &
887	( ztn .LE. sec%ztem(jclass+1))) .OR. &
888	( sec%ztem(jclass) .EQ.99.)) .AND. &
889
890	((( zdep .GE. sec%zlay(jclass)) .AND. &
891	( zdep .LE. sec%zlay(jclass+1))) .OR. &
892	( sec%zlay(jclass) .EQ. 99. )) &
893	)) THEN
894
895	!SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS
896	!----------------------------------------------------------------------------
897	IF (transports_3d(1,jsec,jseg,jk) .GE. 0.0) THEN
898	sec%transport(1,jclass) = sec%transport(1,jclass)+transports_3d(1,jsec,jseg,jk)*1.E-6
899	ELSE
900	sec%transport(2,jclass) = sec%transport(2,jclass)+transports_3d(1,jsec,jseg,jk)*1.E-6
901	ENDIF
902	IF( sec%llstrpond )THEN
903
904	IF ( transports_3d(2,jsec,jseg,jk) .GE. 0.0 ) THEN
905	sec%transport(3,jclass) = sec%transport(3,jclass)+transports_3d(2,jsec,jseg,jk)
906	ELSE
907	sec%transport(4,jclass) = sec%transport(4,jclass)+transports_3d(2,jsec,jseg,jk)
908	ENDIF
909
910	IF ( transports_3d(3,jsec,jseg,jk) .GE. 0.0 ) THEN
911	sec%transport(5,jclass) = sec%transport(5,jclass)+transports_3d(3,jsec,jseg,jk)
912	ELSE
913	sec%transport(6,jclass) = sec%transport(6,jclass)+transports_3d(3,jsec,jseg,jk)
914	ENDIF
915
916	ELSE
917	sec%transport( 3,jclass) = 0._wp
918	sec%transport( 4,jclass) = 0._wp
919	sec%transport( 5,jclass) = 0._wp
920	sec%transport( 6,jclass) = 0._wp
921	ENDIF
922
923	ENDIF ! end of test if point is in class
924
925	END DO ! end of loop on the classes
926
927	END DO ! loop over jk
928
929	#if defined key_si3
930
931	!ICE CASE
932	IF( sec%ll_ice_section )THEN
933
934	IF ( transports_2d(1,jsec,jseg) .GE. 0.0 ) THEN
935	sec%transport( 7,1) = sec%transport( 7,1)+transports_2d(1,jsec,jseg)*1.E-6
936	ELSE
937	sec%transport( 8,1) = sec%transport( 8,1)+transports_2d(1,jsec,jseg)*1.E-6
938	ENDIF
939
940	IF ( transports_2d(3,jsec,jseg) .GE. 0.0 ) THEN
941	sec%transport( 9,1) = sec%transport( 9,1)+transports_2d(2,jsec,jseg)*1.E-6
942	ELSE
943	sec%transport(10,1) = sec%transport(10,1)+transports_2d(2,jsec,jseg)*1.E-6
944	ENDIF
945
946	ENDIF !end of ice case
947	#endif
948
949	END DO !end of loop on the segment
950
951	ELSE !if sec%nb_point =0
952	sec%transport(1:2,:)=0.
953	IF (sec%llstrpond) sec%transport(3:6,:)=0.
954	IF (sec%ll_ice_section) sec%transport(7:10,:)=0.
955	ENDIF !end of sec%nb_point =0 case
956
957	END SUBROUTINE dia_dct_sum
958
959
960	SUBROUTINE dia_dct_wri(kt,ksec,sec)
961	!!-------------------------------------------------------------
962	!! Write transport output in numdct
963	!!
964	!! Purpose: Write transports in ascii files
965	!!
966	!! Method:
967	!! 1. Write volume transports in "volume_transport"
968	!! Unit: Sv : area * Velocity / 1.e6
969	!!
970	!! 2. Write heat transports in "heat_transport"
971	!! Unit: Peta W : area * Velocity * T * rhop * Cp * 1.e-15
972	!!
973	!! 3. Write salt transports in "salt_transport"
974	!! Unit: 10^9 Kg/m^2/s : area * Velocity * S * rhop * 1.e-9
975	!!
976	!!-------------------------------------------------------------
977	!!arguments
978	INTEGER, INTENT(IN) :: kt ! time-step
979	TYPE(SECTION), INTENT(INOUT) :: sec ! section to write
980	INTEGER ,INTENT(IN) :: ksec ! section number
981
982	!!local declarations
983	INTEGER :: jclass ! Dummy loop
984	CHARACTER(len=2) :: classe ! Classname
985	REAL(wp) :: zbnd1,zbnd2 ! Class bounds
986	REAL(wp) :: zslope ! section's slope coeff
987	!
988	REAL(wp), DIMENSION(nb_type_class):: zsumclasses ! 1D workspace
989	!!-------------------------------------------------------------
990
991	zsumclasses(:)=0._wp
992	zslope = sec%slopeSection
993
994
995	DO jclass=1,MAX(1,sec%nb_class-1)
996
997	classe = 'N '
998	zbnd1 = 0._wp
999	zbnd2 = 0._wp
1000	zsumclasses(1:nb_type_class)=zsumclasses(1:nb_type_class)+sec%transport(1:nb_type_class,jclass)
1001
1002
1003	!insitu density classes transports
1004	IF( ( sec%zsigi(jclass) .NE. 99._wp ) .AND. &
1005	( sec%zsigi(jclass+1) .NE. 99._wp ) )THEN
1006	classe = 'DI '
1007	zbnd1 = sec%zsigi(jclass)
1008	zbnd2 = sec%zsigi(jclass+1)
1009	ENDIF
1010	!potential density classes transports
1011	IF( ( sec%zsigp(jclass) .NE. 99._wp ) .AND. &
1012	( sec%zsigp(jclass+1) .NE. 99._wp ) )THEN
1013	classe = 'DP '
1014	zbnd1 = sec%zsigp(jclass)
1015	zbnd2 = sec%zsigp(jclass+1)
1016	ENDIF
1017	!depth classes transports
1018	IF( ( sec%zlay(jclass) .NE. 99._wp ) .AND. &
1019	( sec%zlay(jclass+1) .NE. 99._wp ) )THEN
1020	classe = 'Z '
1021	zbnd1 = sec%zlay(jclass)
1022	zbnd2 = sec%zlay(jclass+1)
1023	ENDIF
1024	!salinity classes transports
1025	IF( ( sec%zsal(jclass) .NE. 99._wp ) .AND. &
1026	( sec%zsal(jclass+1) .NE. 99._wp ) )THEN
1027	classe = 'S '
1028	zbnd1 = sec%zsal(jclass)
1029	zbnd2 = sec%zsal(jclass+1)
1030	ENDIF
1031	!temperature classes transports
1032	IF( ( sec%ztem(jclass) .NE. 99._wp ) .AND. &
1033	( sec%ztem(jclass+1) .NE. 99._wp ) ) THEN
1034	classe = 'T '
1035	zbnd1 = sec%ztem(jclass)
1036	zbnd2 = sec%ztem(jclass+1)
1037	ENDIF
1038
1039	!write volume transport per class
1040	WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, &
1041	jclass,classe,zbnd1,zbnd2,&
1042	sec%transport(1,jclass),sec%transport(2,jclass), &
1043	sec%transport(1,jclass)+sec%transport(2,jclass)
1044
1045	IF( sec%llstrpond )THEN
1046
1047	!write heat transport per class:
1048	WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope, &
1049	jclass,classe,zbnd1,zbnd2,&
1050	sec%transport(3,jclass)1.e-15,sec%transport(4,jclass)1.e-15, &
1051	( sec%transport(3,jclass)+sec%transport(4,jclass) )*1.e-15
1052	!write salt transport per class
1053	WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope, &
1054	jclass,classe,zbnd1,zbnd2,&
1055	sec%transport(5,jclass)1.e-9,sec%transport(6,jclass)1.e-9,&
1056	(sec%transport(5,jclass)+sec%transport(6,jclass))*1.e-9
1057	ENDIF
1058
1059	ENDDO
1060
1061	zbnd1 = 0._wp
1062	zbnd2 = 0._wp
1063	jclass=0
1064
1065	!write total volume transport
1066	WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, &
1067	jclass,"total",zbnd1,zbnd2,&
1068	zsumclasses(1),zsumclasses(2),zsumclasses(1)+zsumclasses(2)
1069
1070	IF( sec%llstrpond )THEN
1071
1072	!write total heat transport
1073	WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope, &
1074	jclass,"total",zbnd1,zbnd2,&
1075	zsumclasses(3)1.e-15,zsumclasses(4)1.e-15,&
1076	(zsumclasses(3)+zsumclasses(4) )*1.e-15
1077	!write total salt transport
1078	WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope, &
1079	jclass,"total",zbnd1,zbnd2,&
1080	zsumclasses(5)1.e-9,zsumclasses(6)1.e-9,&
1081	(zsumclasses(5)+zsumclasses(6))*1.e-9
1082	ENDIF
1083
1084
1085	IF ( sec%ll_ice_section) THEN
1086	!write total ice volume transport
1087	WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,&
1088	jclass,"ice_vol",zbnd1,zbnd2,&
1089	sec%transport(7,1),sec%transport(8,1),&
1090	sec%transport(7,1)+sec%transport(8,1)
1091	!write total ice surface transport
1092	WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,&
1093	jclass,"ice_surf",zbnd1,zbnd2,&
1094	sec%transport(9,1),sec%transport(10,1), &
1095	sec%transport(9,1)+sec%transport(10,1)
1096	ENDIF
1097
1098	118 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3F12.4)
1099	119 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6)
1100	!
1101	END SUBROUTINE dia_dct_wri
1102
1103
1104	FUNCTION interp(ki, kj, kk, cd_point, ptab)
1105	!!----------------------------------------------------------------------
1106	!!
1107	!! Purpose: compute temperature/salinity/density at U-point or V-point
1108	!! --------
1109	!!
1110	!! Method:
1111	!! ------
1112	!!
1113	!! ====> full step and partial step
1114	!!
1115	!!
1116	!! \| I \| I+1 \| Z=temperature/salinity/density at U-poinT
1117	!! \| \| \|
1118	!! ---------------------------------------- 1. Veritcal interpolation: compute zbis
1119	!! \| \| \| interpolation between ptab(I,J,K) and ptab(I,J,K+1)
1120	!! \| \| \| zbis =
1121	!! \| \| \| [ e3w(I+1,J,K)ptab(I,J,K) + ( e3w(I,J,K) - e3w(I+1,J,K) ) ptab(I,J,K-1) ]
1122	!! \| \| \| /[ e3w(I+1,J,K) + e3w(I,J,K) - e3w(I+1,J,K) ]
1123	!! \| \| \|
1124	!! \| \| \| 2. Horizontal interpolation: compute value at U/V point
1125	!!K-1 \| ptab(I,J,K-1) \| \| interpolation between zbis and ptab(I+1,J,K)
1126	!! \| . \| \|
1127	!! \| . \| \| interp = ( 0.5zet2zbis + 0.5zet1ptab(I+1,J,K) )/(0.5zet2+0.5zet1)
1128	!! \| . \| \|
1129	!! ------------------------------------------
1130	!! \| . \| \|
1131	!! \| . \| \|
1132	!! \| . \| \|
1133	!!K \| zbis.......U...ptab(I+1,J,K) \|
1134	!! \| . \| \|
1135	!! \| ptab(I,J,K) \| \|
1136	!! \| \|------------------\|
1137	!! \| \| partials \|
1138	!! \| \| steps \|
1139	!! -------------------------------------------
1140	!! <----zet1------><----zet2--------->
1141	!!
1142	!!
1143	!! ====> s-coordinate
1144	!!
1145	!! \| \| \| 1. Compute distance between T1 and U points: SQRT( zdep1^2 + (0.5 * zet1 )^2
1146	!! \| \| \| Compute distance between T2 and U points: SQRT( zdep2^2 + (0.5 * zet2 )^2
1147	!! \| \| ptab(I+1,J,K) \|
1148	!! \| \| T2 \| 2. Interpolation between T1 and T2 values at U point
1149	!! \| \| ^ \|
1150	!! \| \| \| zdep2 \|
1151	!! \| \| \| \|
1152	!! \| ^ U v \|
1153	!! \| \| \| \|
1154	!! \| \| zdep1 \| \|
1155	!! \| v \| \|
1156	!! \| T1 \| \|
1157	!! \| ptab(I,J,K) \| \|
1158	!! \| \| \|
1159	!! \| \| \|
1160	!!
1161	!! <----zet1--------><----zet2--------->
1162	!!
1163	!!----------------------------------------------------------------------
1164	!*arguments
1165	INTEGER, INTENT(IN) :: ki, kj, kk ! coordinate of point
1166	CHARACTER(len=1), INTENT(IN) :: cd_point ! type of point (U, V)
1167	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ptab ! variable to compute at (ki, kj, kk )
1168	REAL(wp) :: interp ! interpolated variable
1169
1170	!*local declations
1171	INTEGER :: ii1, ij1, ii2, ij2 ! local integer
1172	REAL(wp):: ze3t, ze3, zwgt1, zwgt2, zbis, zdepu ! local real
1173	REAL(wp):: zet1, zet2 ! weight for interpolation
1174	REAL(wp):: zdep1,zdep2 ! differences of depth
1175	REAL(wp):: zmsk ! mask value
1176	!!----------------------------------------------------------------------
1177
1178	IF( cd_point=='U' )THEN
1179	ii1 = ki ; ij1 = kj
1180	ii2 = ki+1 ; ij2 = kj
1181
1182	zet1=e1t(ii1,ij1)
1183	zet2=e1t(ii2,ij2)
1184	zmsk=umask(ii1,ij1,kk)
1185
1186
1187	ELSE ! cd_point=='V'
1188	ii1 = ki ; ij1 = kj
1189	ii2 = ki ; ij2 = kj+1
1190
1191	zet1=e2t(ii1,ij1)
1192	zet2=e2t(ii2,ij2)
1193	zmsk=vmask(ii1,ij1,kk)
1194
1195	ENDIF
1196
1197	IF( ln_sco )THEN ! s-coordinate case
1198
1199	zdepu = ( gdept_n(ii1,ij1,kk) + gdept_n(ii2,ij2,kk) ) * 0.5_wp
1200	zdep1 = gdept_n(ii1,ij1,kk) - zdepu
1201	zdep2 = gdept_n(ii2,ij2,kk) - zdepu
1202
1203	! weights
1204	zwgt1 = SQRT( ( 0.5 * zet1 ) * ( 0.5 * zet1 ) + ( zdep1 * zdep1 ) )
1205	zwgt2 = SQRT( ( 0.5 * zet2 ) * ( 0.5 * zet2 ) + ( zdep2 * zdep2 ) )
1206
1207	! result
1208	interp = zmsk * ( zwgt2 * ptab(ii1,ij1,kk) + zwgt1 * ptab(ii1,ij1,kk) ) / ( zwgt2 + zwgt1 )
1209
1210
1211	ELSE ! full step or partial step case
1212
1213	ze3t = e3t_n(ii2,ij2,kk) - e3t_n(ii1,ij1,kk)
1214	zwgt1 = ( e3w_n(ii2,ij2,kk) - e3w_n(ii1,ij1,kk) ) / e3w_n(ii2,ij2,kk)
1215	zwgt2 = ( e3w_n(ii1,ij1,kk) - e3w_n(ii2,ij2,kk) ) / e3w_n(ii1,ij1,kk)
1216
1217	IF(kk .NE. 1)THEN
1218
1219	IF( ze3t >= 0. )THEN
1220	! zbis
1221	zbis = ptab(ii2,ij2,kk) + zwgt1 * ( ptab(ii2,ij2,kk-1) - ptab(ii2,ij2,kk) )
1222	! result
1223	interp = zmsk * ( zet2 * ptab(ii1,ij1,kk) + zet1 * zbis )/( zet1 + zet2 )
1224	ELSE
1225	! zbis
1226	zbis = ptab(ii1,ij1,kk) + zwgt2 * ( ptab(ii1,ij1,kk-1) - ptab(ii1,ij2,kk) )
1227	! result
1228	interp = zmsk * ( zet2 * zbis + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 )
1229	ENDIF
1230
1231	ELSE
1232	interp = zmsk * ( zet2 * ptab(ii1,ij1,kk) + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 )
1233	ENDIF
1234
1235	ENDIF
1236	!
1237	END FUNCTION interp
1238
1239	#else
1240	!!----------------------------------------------------------------------
1241	!! Default option : Dummy module
1242	!!----------------------------------------------------------------------
1243	LOGICAL, PUBLIC, PARAMETER :: lk_diadct = .FALSE. !: diamht flag
1244	PUBLIC
1245	!! $Id$
1246	CONTAINS
1247
1248	SUBROUTINE dia_dct_init ! Dummy routine
1249	WRITE(,) 'dia_dct_init: You should not have seen this print! error?', kt
1250	END SUBROUTINE dia_dct_init
1251
1252	SUBROUTINE dia_dct( kt ) ! Dummy routine
1253	INTEGER, INTENT( in ) :: kt ! ocean time-step index
1254	WRITE(,) 'dia_dct: You should not have seen this print! error?', kt
1255	END SUBROUTINE dia_dct
1256	#endif
1257
1258	!!======================================================================
1259	END MODULE diadct

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

Download in other formats: