Context Navigation

diadct.F90 @ 10986

Last change on this file since 10986 was 10986, checked in by andmirek, 5 years ago
GMED 462 add flush
File size: 59.0 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 license (see ./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_stop( 'STOP', '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 .AND. nprint > 2) 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	IF(lflush) CALL FLUSH(numout)
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	IF(lflush) CALL FLUSH(numout)
225	ENDIF
226
227
228	! Compute transport and write only at nn_dctwri
229	IF( MOD(kt,nn_dct)==0 ) THEN
230
231	DO jsec=1,nb_sec
232
233	!debug this section computing ?
234	lldebug=.FALSE.
235	IF( (jsec==nn_secdebug .OR. nn_secdebug==-1) .AND. kt==nit000+nn_dct-1 ) lldebug=.TRUE.
236
237	!Compute transport through section
238	CALL transport(secs(jsec),lldebug,jsec)
239
240	ENDDO
241
242	IF( MOD(kt,nn_dctwri)==0 )THEN
243
244	IF( kt==nit000+nn_dctwri-1 ) THEN
245	WRITE(numout,*)" diadct: average transports and write at kt = ",kt
246	IF(lflush) CALL FLUSH(numout)
247	ENDIF
248
249	!! divide arrays by nn_dctwri/nn_dct to obtain average
250	transports_3d(:,:,:,:)=transports_3d(:,:,:,:)/(nn_dctwri/nn_dct)
251	transports_2d(:,:,:) =transports_2d(:,:,:) /(nn_dctwri/nn_dct)
252
253	! Sum over each class
254	DO jsec=1,nb_sec
255	CALL dia_dct_sum(secs(jsec),jsec)
256	ENDDO
257
258	!Sum on all procs
259	IF( lk_mpp )THEN
260	ish(1) = nb_sec_maxnb_type_classnb_class_max
261	ish2 = (/nb_sec_max,nb_type_class,nb_class_max/)
262	DO jsec=1,nb_sec ; zsum(jsec,:,:) = secs(jsec)%transport(:,:) ; ENDDO
263	zwork(:)= RESHAPE(zsum(:,:,:), ish )
264	CALL mpp_sum('diadct', zwork, ish(1))
265	zsum(:,:,:)= RESHAPE(zwork,ish2)
266	DO jsec=1,nb_sec ; secs(jsec)%transport(:,:) = zsum(jsec,:,:) ; ENDDO
267	ENDIF
268
269	!Write the transport
270	DO jsec=1,nb_sec
271
272	IF( lwm )CALL dia_dct_wri(kt,jsec,secs(jsec))
273
274	!nullify transports values after writing
275	transports_3d(:,jsec,:,:)=0.
276	transports_2d(:,jsec,: )=0.
277	secs(jsec)%transport(:,:)=0.
278
279	ENDDO
280
281	ENDIF
282
283	ENDIF
284
285	IF( lk_mpp )THEN
286	itotal = nb_sec_maxnb_type_classnb_class_max
287	DEALLOCATE( zwork , zsum )
288	ENDIF
289
290	IF( ln_timing ) CALL timing_stop('dia_dct')
291	!
292	END SUBROUTINE dia_dct
293
294
295	SUBROUTINE readsec
296	!!---------------------------------------------------------------------
297	!! * ROUTINE readsec *
298	!!
299	!! ** Purpose:
300	!! Read a binary file(section_ijglobal.diadct)
301	!! generated by the tools "NEMOGCM/TOOLS/SECTIONS_DIADCT"
302	!!
303	!!
304	!!---------------------------------------------------------------------
305	INTEGER :: iptglo , iptloc ! Global and local number of points for a section
306	INTEGER :: isec, iiglo, ijglo, iiloc, ijloc,iost,i1 ,i2 ! temporary integer
307	INTEGER :: jsec, jpt ! dummy loop indices
308	INTEGER, DIMENSION(2) :: icoord
309	LOGICAL :: llbon, lldebug ! local logical
310	CHARACTER(len=160) :: clname ! filename
311	CHARACTER(len=200) :: cltmp
312	CHARACTER(len=200) :: clformat !automatic format
313	TYPE(POINT_SECTION),DIMENSION(nb_point_max) ::coordtemp !contains listpoints coordinates read in the file
314	INTEGER, DIMENSION(nb_point_max) :: directemp !contains listpoints directions read in the files
315	!!-------------------------------------------------------------------------------------
316
317	!open input file
318	!---------------
319	CALL ctl_opn( numdct_in, 'section_ijglobal.diadct', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
320
321	!---------------
322	!Read input file
323	!---------------
324
325	DO jsec=1,nb_sec_max !loop on the nb_sec sections
326
327	IF ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) THEN
328	WRITE(numout,*)'debuging for section number: ',jsec
329	IF(lflush) CALL FLUSH(numout)
330	ENDIF
331
332	!initialization
333	!---------------
334	secs(jsec)%name=''
335	secs(jsec)%llstrpond = .FALSE. ; secs(jsec)%ll_ice_section = .FALSE.
336	secs(jsec)%ll_date_line = .FALSE. ; secs(jsec)%nb_class = 0
337	secs(jsec)%zsigi = 99._wp ; secs(jsec)%zsigp = 99._wp
338	secs(jsec)%zsal = 99._wp ; secs(jsec)%ztem = 99._wp
339	secs(jsec)%zlay = 99._wp
340	secs(jsec)%transport = 0._wp ; secs(jsec)%nb_point = 0
341
342	!read section's number / name / computing choices / classes / slopeSection / points number
343	!-----------------------------------------------------------------------------------------
344	READ(numdct_in,iostat=iost)isec
345	IF (iost .NE. 0 )EXIT !end of file
346	WRITE(cltmp,'(a,i4.4,a,i4.4)')'diadct: read sections : Problem of section number: isec= ',isec,' and jsec= ',jsec
347	IF( jsec .NE. isec ) CALL ctl_stop( cltmp )
348
349	IF( jsec==nn_secdebug .OR. nn_secdebug==-1 ) THEN
350	WRITE(numout,*)"isec ",isec
351	IF(lflush) CALL FLUSH(numout)
352	ENDIF
353
354	READ(numdct_in)secs(jsec)%name
355	READ(numdct_in)secs(jsec)%llstrpond
356	READ(numdct_in)secs(jsec)%ll_ice_section
357	READ(numdct_in)secs(jsec)%ll_date_line
358	READ(numdct_in)secs(jsec)%coordSec
359	READ(numdct_in)secs(jsec)%nb_class
360	READ(numdct_in)secs(jsec)%zsigi
361	READ(numdct_in)secs(jsec)%zsigp
362	READ(numdct_in)secs(jsec)%zsal
363	READ(numdct_in)secs(jsec)%ztem
364	READ(numdct_in)secs(jsec)%zlay
365	READ(numdct_in)secs(jsec)%slopeSection
366	READ(numdct_in)iptglo
367
368	!debug
369	!-----
370
371	IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )THEN
372
373	WRITE(clformat,'(a,i2,a)') '(A40,', nb_class_max,'(f8.3,1X))'
374
375	WRITE(numout,*) " Section name : ",TRIM(secs(jsec)%name)
376	WRITE(numout,*) " Compute heat and salt transport ? ",secs(jsec)%llstrpond
377	WRITE(numout,*) " Compute ice transport ? ",secs(jsec)%ll_ice_section
378	WRITE(numout,*) " Section crosses date-line ? ",secs(jsec)%ll_date_line
379	WRITE(numout,*) " Slope section : ",secs(jsec)%slopeSection
380	WRITE(numout,*) " Number of points in the section: ",iptglo
381	WRITE(numout,*) " Number of classes ",secs(jsec)%nb_class
382	WRITE(numout,clformat)" Insitu density classes : ",secs(jsec)%zsigi
383	WRITE(numout,clformat)" Potential density classes : ",secs(jsec)%zsigp
384	WRITE(numout,clformat)" Salinity classes : ",secs(jsec)%zsal
385	WRITE(numout,clformat)" Temperature classes : ",secs(jsec)%ztem
386	WRITE(numout,clformat)" Depth classes : ",secs(jsec)%zlay
387	IF(lflush) CALL FLUSH(numout)
388	ENDIF
389
390	IF( iptglo /= 0 )THEN
391
392	!read points'coordinates and directions
393	!--------------------------------------
394	coordtemp(:) = POINT_SECTION(0,0) !list of points read
395	directemp(:) = 0 !value of directions of each points
396	DO jpt=1,iptglo
397	READ(numdct_in) i1, i2
398	coordtemp(jpt)%I = i1
399	coordtemp(jpt)%J = i2
400	ENDDO
401	READ(numdct_in) directemp(1:iptglo)
402
403	!debug
404	!-----
405	IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )THEN
406	WRITE(numout,*)" List of points in global domain:"
407	DO jpt=1,iptglo
408	WRITE(numout,*)' # I J ',jpt,coordtemp(jpt),directemp(jpt)
409	ENDDO
410	IF(lflush) CALL FLUSH(numout)
411	ENDIF
412
413	!Now each proc selects only points that are in its domain:
414	!--------------------------------------------------------
415	iptloc = 0 ! initialize number of points selected
416	DO jpt = 1, iptglo ! loop on listpoint read in the file
417	!
418	iiglo=coordtemp(jpt)%I ! global coordinates of the point
419	ijglo=coordtemp(jpt)%J ! "
420
421	IF( iiglo==jpiglo .AND. nimpp==1 ) iiglo = 2 !!gm BUG: Hard coded periodicity !
422
423	iiloc=iiglo-nimpp+1 ! local coordinates of the point
424	ijloc=ijglo-njmpp+1 ! "
425
426	!verify if the point is on the local domain:(1,nlei)*(1,nlej)
427	IF( iiloc >= 1 .AND. iiloc <= nlei .AND. &
428	ijloc >= 1 .AND. ijloc <= nlej )THEN
429	iptloc = iptloc + 1 ! count local points
430	secs(jsec)%listPoint(iptloc) = POINT_SECTION(mi0(iiglo),mj0(ijglo)) ! store local coordinates
431	secs(jsec)%direction(iptloc) = directemp(jpt) ! store local direction
432	ENDIF
433	!
434	END DO
435
436	secs(jsec)%nb_point=iptloc !store number of section's points
437
438	!debug
439	!-----
440	IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )THEN
441	WRITE(numout,*)" List of points selected by the proc:"
442	DO jpt = 1,iptloc
443	iiglo = secs(jsec)%listPoint(jpt)%I + nimpp - 1
444	ijglo = secs(jsec)%listPoint(jpt)%J + njmpp - 1
445	WRITE(numout,*)' # I J : ',iiglo,ijglo
446	ENDDO
447	IF(lflush) CALL FLUSH(numout)
448	ENDIF
449
450	IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN
451	DO jpt = 1,iptloc
452	iiglo = secs(jsec)%listPoint(jpt)%I + nimpp - 1
453	ijglo = secs(jsec)%listPoint(jpt)%J + njmpp - 1
454	ENDDO
455	ENDIF
456
457	!remove redundant points between processors
458	!------------------------------------------
459	lldebug = .FALSE. ; IF ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) lldebug = .TRUE.
460	IF( iptloc .NE. 0 )THEN
461	CALL removepoints(secs(jsec),'I','top_list',lldebug)
462	CALL removepoints(secs(jsec),'I','bot_list',lldebug)
463	CALL removepoints(secs(jsec),'J','top_list',lldebug)
464	CALL removepoints(secs(jsec),'J','bot_list',lldebug)
465	ENDIF
466	IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN
467	DO jpt = 1,secs(jsec)%nb_point
468	iiglo = secs(jsec)%listPoint(jpt)%I + nimpp - 1
469	ijglo = secs(jsec)%listPoint(jpt)%J + njmpp - 1
470	ENDDO
471	ENDIF
472
473	!debug
474	!-----
475	IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )THEN
476	WRITE(numout,*)" List of points after removepoints:"
477	iptloc = secs(jsec)%nb_point
478	DO jpt = 1,iptloc
479	iiglo = secs(jsec)%listPoint(jpt)%I + nimpp - 1
480	ijglo = secs(jsec)%listPoint(jpt)%J + njmpp - 1
481	WRITE(numout,*)' # I J : ',iiglo,ijglo
482	IF(lflush) CALL FLUSH(numout)
483	ENDDO
484	ENDIF
485
486	ELSE ! iptglo = 0
487	IF( jsec==nn_secdebug .OR. nn_secdebug==-1 ) THEN
488	WRITE(numout,*)' No points for this section.'
489	IF(lflush) CALL FLUSH(numout)
490	ENDIF
491	ENDIF
492
493	ENDDO !end of the loop on jsec
494
495	nb_sec = jsec-1 !number of section read in the file
496	!
497	END SUBROUTINE readsec
498
499
500	SUBROUTINE removepoints(sec,cdind,cdextr,ld_debug)
501	!!---------------------------------------------------------------------------
502	!! *** function removepoints
503	!!
504	!! ** Purpose :: Remove points which are common to 2 procs
505	!!
506	!----------------------------------------------------------------------------
507	!! * arguments
508	TYPE(SECTION),INTENT(INOUT) :: sec
509	CHARACTER(len=1),INTENT(IN) :: cdind ! = 'I'/'J'
510	CHARACTER(len=8),INTENT(IN) :: cdextr ! = 'top_list'/'bot_list'
511	LOGICAL,INTENT(IN) :: ld_debug
512
513	!! * Local variables
514	INTEGER :: iextr ,& !extremity of listpoint that we verify
515	iind ,& !coord of listpoint that we verify
516	itest ,& !indice value of the side of the domain
517	!where points could be redundant
518	isgn ,& ! isgn= 1 : scan listpoint from start to end
519	! isgn=-1 : scan listpoint from end to start
520	istart,iend !first and last points selected in listpoint
521	INTEGER :: jpoint !loop on list points
522	INTEGER, DIMENSION(nb_point_max) :: idirec !contains temporary sec%direction
523	INTEGER, DIMENSION(2,nb_point_max) :: icoord !contains temporary sec%listpoint
524	!----------------------------------------------------------------------------
525	!
526	IF( ld_debug ) THEN
527	WRITE(numout,*)' -------------------------'
528	WRITE(numout,*)' removepoints in listpoint'
529	IF(lflush) CALL FLUSH(numout)
530	ENDIF
531
532	!iextr=extremity of list_point that we verify
533	IF ( cdextr=='bot_list' )THEN ; iextr=1 ; isgn=1
534	ELSE IF ( cdextr=='top_list' )THEN ; iextr=sec%nb_point ; isgn=-1
535	ELSE ; CALL ctl_stop("removepoints :Wrong value for cdextr")
536	ENDIF
537
538	!which coordinate shall we verify ?
539	IF ( cdind=='I' )THEN ; itest=nlei ; iind=1
540	ELSE IF ( cdind=='J' )THEN ; itest=nlej ; iind=2
541	ELSE ; CALL ctl_stop("removepoints :Wrong value for cdind")
542	ENDIF
543
544	IF( ld_debug )THEN
545	WRITE(numout,*)' case: coord/list extr/domain side'
546	WRITE(numout,*)' ', cdind,' ',cdextr,' ',itest
547	WRITE(numout,*)' Actual number of points: ',sec%nb_point
548	IF(lflush) CALL FLUSH(numout)
549	ENDIF
550
551	icoord(1,1:nb_point_max) = sec%listPoint%I
552	icoord(2,1:nb_point_max) = sec%listPoint%J
553	idirec = sec%direction
554	sec%listPoint = POINT_SECTION(0,0)
555	sec%direction = 0
556
557	jpoint=iextr+isgn
558	DO WHILE( jpoint .GE. 1 .AND. jpoint .LE. sec%nb_point )
559	IF( icoord( iind,jpoint-isgn ) == itest .AND. icoord( iind,jpoint ) == itest )THEN ; jpoint=jpoint+isgn
560	ELSE ; EXIT
561	ENDIF
562	ENDDO
563
564	IF( cdextr=='bot_list')THEN ; istart=jpoint-1 ; iend=sec%nb_point
565	ELSE ; istart=1 ; iend=jpoint+1
566	ENDIF
567
568	sec%listPoint(1:1+iend-istart)%I = icoord(1,istart:iend)
569	sec%listPoint(1:1+iend-istart)%J = icoord(2,istart:iend)
570	sec%direction(1:1+iend-istart) = idirec(istart:iend)
571	sec%nb_point = iend-istart+1
572
573	IF( ld_debug )THEN
574	WRITE(numout,*)' Number of points after removepoints :',sec%nb_point
575	WRITE(numout,*)' sec%direction after removepoints :',sec%direction(1:sec%nb_point)
576	IF(lflush) CALL FLUSH(numout)
577	ENDIF
578	!
579	END SUBROUTINE removepoints
580
581
582	SUBROUTINE transport(sec,ld_debug,jsec)
583	!!-------------------------------------------------------------------------------------------
584	!! * ROUTINE transport *
585	!!
586	!! Purpose :: Compute the transport for each point in a section
587	!!
588	!! Method :: Loop over each segment, and each vertical level and add the transport
589	!! Be aware :
590	!! One section is a sum of segments
591	!! One segment is defined by 2 consecutive points in sec%listPoint
592	!! All points of sec%listPoint are positioned on the F-point of the cell
593	!!
594	!! There are two loops:
595	!! loop on the segment between 2 nodes
596	!! loop on the level jk !!
597	!!
598	!! Output :: Arrays containing the volume,density,heat,salt transports for each i
599	!! point in a section, summed over each nn_dct.
600	!!
601	!!-------------------------------------------------------------------------------------------
602	TYPE(SECTION),INTENT(INOUT) :: sec
603	LOGICAL ,INTENT(IN) :: ld_debug
604	INTEGER ,INTENT(IN) :: jsec ! numeric identifier of section
605	!
606	INTEGER :: jk, jseg, jclass,jl, isgnu, isgnv ! loop on level/segment/classes/ice categories
607	REAL(wp):: zumid, zvmid, zumid_ice, zvmid_ice ! U/V ocean & ice velocity on a cell segment
608	REAL(wp):: zTnorm ! transport of velocity through one cell's sides
609	REAL(wp):: ztn, zsn, zrhoi, zrhop, zsshn, zdep ! temperature/salinity/potential density/ssh/depth at u/v point
610	TYPE(POINT_SECTION) :: k
611	!!--------------------------------------------------------
612	!
613	IF( ld_debug ) THEN
614	WRITE(numout,*)' Compute transport'
615	IF(lflush) CALL FLUSH(numout)
616	ENDIF
617
618	!---------------------------!
619	! COMPUTE TRANSPORT !
620	!---------------------------!
621	IF(sec%nb_point .NE. 0)THEN
622
623	!----------------------------------------------------------------------------------------------------
624	!Compute sign for velocities:
625	!
626	!convention:
627	! non horizontal section: direction + is toward left hand of section
628	! horizontal section: direction + is toward north of section
629	!
630	!
631	! slopeSection < 0 slopeSection > 0 slopeSection=inf slopeSection=0
632	! ---------------- ----------------- --------------- --------------
633	!
634	! isgnv=1 direction +
635	! ______ _____ ______
636	! \| //\| \| \| direction +
637	! \| isgnu=1 // \| \|isgnu=1 \|isgnu=1 /\|\
638	! \|_______ // ______\| \\ \| ---\ \|
639	! \| \| isgnv=-1 \\ \| \| ---/ direction + ____________
640	! \| \| __\\\| \|
641	! \| \| direction + \| isgnv=1
642	!
643	!----------------------------------------------------------------------------------------------------
644	isgnu = 1
645	IF( sec%slopeSection .GT. 0 ) THEN ; isgnv = -1
646	ELSE ; isgnv = 1
647	ENDIF
648	IF( sec%slopeSection .GE. 9999. ) isgnv = 1
649
650	IF( ld_debug ) THEN
651	write(numout,*)"sec%slopeSection isgnu isgnv ",sec%slopeSection,isgnu,isgnv
652	IF(lflush) CALL FLUSH(numout)
653	ENDIF
654
655	!--------------------------------------!
656	! LOOP ON THE SEGMENT BETWEEN 2 NODES !
657	!--------------------------------------!
658	DO jseg=1,MAX(sec%nb_point-1,0)
659
660	!-------------------------------------------------------------------------------------------
661	! Select the appropriate coordinate for computing the velocity of the segment
662	!
663	! CASE(0) Case (2)
664	! ------- --------
665	! listPoint(jseg) listPoint(jseg+1) listPoint(jseg) F(i,j)
666	! F(i,j)----------V(i+1,j)-------F(i+1,j) \|
667	! \|
668	! \|
669	! \|
670	! Case (3) U(i,j)
671	! -------- \|
672	! \|
673	! listPoint(jseg+1) F(i,j+1) \|
674	! \| \|
675	! \| \|
676	! \| listPoint(jseg+1) F(i,j-1)
677	! \|
678	! \|
679	! U(i,j+1)
680	! \| Case(1)
681	! \| ------
682	! \|
683	! \| listPoint(jseg+1) listPoint(jseg)
684	! \| F(i-1,j)-----------V(i,j) -------f(jseg)
685	! listPoint(jseg) F(i,j)
686	!
687	!-------------------------------------------------------------------------------------------
688
689	SELECT CASE( sec%direction(jseg) )
690	CASE(0) ; k = sec%listPoint(jseg)
691	CASE(1) ; k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
692	CASE(2) ; k = sec%listPoint(jseg)
693	CASE(3) ; k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
694	END SELECT
695
696	!---------------------------\|
697	! LOOP ON THE LEVEL \|
698	!---------------------------\|
699	DO jk = 1, mbkt(k%I,k%J) !Sum of the transport on the vertical
700	! ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
701	SELECT CASE( sec%direction(jseg) )
702	CASE(0,1)
703	ztn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) )
704	zsn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) )
705	zrhop = interp(k%I,k%J,jk,'V',rhop)
706	zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0)
707	zsshn = 0.5( sshn(k%I,k%J) + sshn(k%I,k%J+1) ) vmask(k%I,k%J,1)
708	CASE(2,3)
709	ztn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) )
710	zsn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) )
711	zrhop = interp(k%I,k%J,jk,'U',rhop)
712	zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0)
713	zsshn = 0.5( sshn(k%I,k%J) + sshn(k%I+1,k%J) ) umask(k%I,k%J,1)
714	END SELECT
715	!
716	zdep= gdept_n(k%I,k%J,jk)
717
718	SELECT CASE( sec%direction(jseg) ) !compute velocity with the correct direction
719	CASE(0,1)
720	zumid=0._wp
721	zvmid=isgnvvn(k%I,k%J,jk)vmask(k%I,k%J,jk)
722	CASE(2,3)
723	zumid=isgnuun(k%I,k%J,jk)umask(k%I,k%J,jk)
724	zvmid=0._wp
725	END SELECT
726
727	!zTnorm=transport through one cell;
728	!velocity* cell's length * cell's thickness
729	zTnorm = zumide2u(k%I,k%J) e3u_n(k%I,k%J,jk) &
730	& + zvmide1v(k%I,k%J) e3v_n(k%I,k%J,jk)
731
732	!!gm THIS is WRONG no transport due to ssh in linear free surface case !!!!!
733	IF( ln_linssh ) THEN !add transport due to free surface
734	IF( jk==1 ) THEN
735	zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) &
736	& + zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk)
737	ENDIF
738	ENDIF
739	!!gm end
740	!COMPUTE TRANSPORT
741
742	transports_3d(1,jsec,jseg,jk) = transports_3d(1,jsec,jseg,jk) + zTnorm
743
744	IF( sec%llstrpond ) THEN
745	transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk) + zTnorm * ztn * zrhop * rcp
746	transports_3d(3,jsec,jseg,jk) = transports_3d(3,jsec,jseg,jk) + zTnorm * zsn * zrhop * 0.001
747	ENDIF
748
749	END DO !end of loop on the level
750
751	#if defined key_si3
752
753	!ICE CASE
754	!------------
755	IF( sec%ll_ice_section )THEN
756	SELECT CASE (sec%direction(jseg))
757	CASE(0)
758	zumid_ice = 0
759	zvmid_ice = isgnv0.5(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
760	CASE(1)
761	zumid_ice = 0
762	zvmid_ice = isgnv0.5(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
763	CASE(2)
764	zvmid_ice = 0
765	zumid_ice = isgnu0.5(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
766	CASE(3)
767	zvmid_ice = 0
768	zumid_ice = isgnu0.5(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
769	END SELECT
770
771	zTnorm=zumid_icee2u(k%I,k%J)+zvmid_icee1v(k%I,k%J)
772
773	#if defined key_si3
774	DO jl=1,jpl
775	transports_2d(1,jsec,jseg) = transports_2d(1,jsec,jseg) + (zTnorm)* &
776	a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) * &
777	( h_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) + &
778	h_s(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) )
779
780	transports_2d(2,jsec,jseg) = transports_2d(2,jsec,jseg) + (zTnorm)* &
781	a_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl)
782	END DO
783	#endif
784
785	ENDIF !end of ice case
786	#endif
787
788	END DO !end of loop on the segment
789
790	ENDIF !end of sec%nb_point =0 case
791	!
792	END SUBROUTINE transport
793
794
795	SUBROUTINE dia_dct_sum(sec,jsec)
796	!!-------------------------------------------------------------
797	!! Purpose: Average the transport over nn_dctwri time steps
798	!! and sum over the density/salinity/temperature/depth classes
799	!!
800	!! Method: Sum over relevant grid cells to obtain values
801	!! for each class
802	!! There are several loops:
803	!! loop on the segment between 2 nodes
804	!! loop on the level jk
805	!! loop on the density/temperature/salinity/level classes
806	!! test on the density/temperature/salinity/level
807	!!
808	!! Note: Transport through a given section is equal to the sum of transports
809	!! computed on each proc.
810	!! On each proc,transport is equal to the sum of transport computed through
811	!! segments linking each point of sec%listPoint with the next one.
812	!!
813	!!-------------------------------------------------------------
814	TYPE(SECTION),INTENT(INOUT) :: sec
815	INTEGER ,INTENT(IN) :: jsec ! numeric identifier of section
816
817	TYPE(POINT_SECTION) :: k
818	INTEGER :: jk,jseg,jclass ! dummy variables for looping on level/segment/classes
819	REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zdep ! temperature/salinity/ssh/potential density /depth at u/v point
820	!!-------------------------------------------------------------
821
822	!! Sum the relevant segments to obtain values for each class
823	IF(sec%nb_point .NE. 0)THEN
824
825	!--------------------------------------!
826	! LOOP ON THE SEGMENT BETWEEN 2 NODES !
827	!--------------------------------------!
828	DO jseg=1,MAX(sec%nb_point-1,0)
829
830	!-------------------------------------------------------------------------------------------
831	! Select the appropriate coordinate for computing the velocity of the segment
832	!
833	! CASE(0) Case (2)
834	! ------- --------
835	! listPoint(jseg) listPoint(jseg+1) listPoint(jseg) F(i,j)
836	! F(i,j)----------V(i+1,j)-------F(i+1,j) \|
837	! \|
838	! \|
839	! \|
840	! Case (3) U(i,j)
841	! -------- \|
842	! \|
843	! listPoint(jseg+1) F(i,j+1) \|
844	! \| \|
845	! \| \|
846	! \| listPoint(jseg+1) F(i,j-1)
847	! \|
848	! \|
849	! U(i,j+1)
850	! \| Case(1)
851	! \| ------
852	! \|
853	! \| listPoint(jseg+1) listPoint(jseg)
854	! \| F(i-1,j)-----------V(i,j) -------f(jseg)
855	! listPoint(jseg) F(i,j)
856	!
857	!-------------------------------------------------------------------------------------------
858
859	SELECT CASE( sec%direction(jseg) )
860	CASE(0) ; k = sec%listPoint(jseg)
861	CASE(1) ; k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
862	CASE(2) ; k = sec%listPoint(jseg)
863	CASE(3) ; k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
864	END SELECT
865
866	!---------------------------\|
867	! LOOP ON THE LEVEL \|
868	!---------------------------\|
869	!Sum of the transport on the vertical
870	DO jk=1,mbkt(k%I,k%J)
871
872	! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
873	SELECT CASE( sec%direction(jseg) )
874	CASE(0,1)
875	ztn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) )
876	zsn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) )
877	zrhop = interp(k%I,k%J,jk,'V',rhop)
878	zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0)
879
880	CASE(2,3)
881	ztn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) )
882	zsn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) )
883	zrhop = interp(k%I,k%J,jk,'U',rhop)
884	zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0)
885	zsshn = 0.5( sshn(k%I,k%J) + sshn(k%I+1,k%J) ) umask(k%I,k%J,1)
886	END SELECT
887
888	zdep= gdept_n(k%I,k%J,jk)
889
890	!-------------------------------
891	! LOOP ON THE DENSITY CLASSES \|
892	!-------------------------------
893	!The computation is made for each density/temperature/salinity/depth class
894	DO jclass=1,MAX(1,sec%nb_class-1)
895
896	!----------------------------------------------!
897	!TEST ON THE DENSITY/SALINITY/TEMPERATURE/LEVEL!
898	!----------------------------------------------!
899
900	IF ( ( &
901	((( zrhop .GE. (sec%zsigp(jclass)+1000. )) .AND. &
902	( zrhop .LE. (sec%zsigp(jclass+1)+1000. ))) .OR. &
903	( sec%zsigp(jclass) .EQ. 99.)) .AND. &
904
905	((( zrhoi .GE. (sec%zsigi(jclass) + 1000. )) .AND. &
906	( zrhoi .LE. (sec%zsigi(jclass+1)+1000. ))) .OR. &
907	( sec%zsigi(jclass) .EQ. 99.)) .AND. &
908
909	((( zsn .GT. sec%zsal(jclass)) .AND. &
910	( zsn .LE. sec%zsal(jclass+1))) .OR. &
911	( sec%zsal(jclass) .EQ. 99.)) .AND. &
912
913	((( ztn .GE. sec%ztem(jclass)) .AND. &
914	( ztn .LE. sec%ztem(jclass+1))) .OR. &
915	( sec%ztem(jclass) .EQ.99.)) .AND. &
916
917	((( zdep .GE. sec%zlay(jclass)) .AND. &
918	( zdep .LE. sec%zlay(jclass+1))) .OR. &
919	( sec%zlay(jclass) .EQ. 99. )) &
920	)) THEN
921
922	!SUM THE TRANSPORTS FOR EACH CLASSES FOR THE POSITIVE AND NEGATIVE DIRECTIONS
923	!----------------------------------------------------------------------------
924	IF (transports_3d(1,jsec,jseg,jk) .GE. 0.0) THEN
925	sec%transport(1,jclass) = sec%transport(1,jclass)+transports_3d(1,jsec,jseg,jk)*1.E-6
926	ELSE
927	sec%transport(2,jclass) = sec%transport(2,jclass)+transports_3d(1,jsec,jseg,jk)*1.E-6
928	ENDIF
929	IF( sec%llstrpond )THEN
930
931	IF ( transports_3d(2,jsec,jseg,jk) .GE. 0.0 ) THEN
932	sec%transport(3,jclass) = sec%transport(3,jclass)+transports_3d(2,jsec,jseg,jk)
933	ELSE
934	sec%transport(4,jclass) = sec%transport(4,jclass)+transports_3d(2,jsec,jseg,jk)
935	ENDIF
936
937	IF ( transports_3d(3,jsec,jseg,jk) .GE. 0.0 ) THEN
938	sec%transport(5,jclass) = sec%transport(5,jclass)+transports_3d(3,jsec,jseg,jk)
939	ELSE
940	sec%transport(6,jclass) = sec%transport(6,jclass)+transports_3d(3,jsec,jseg,jk)
941	ENDIF
942
943	ELSE
944	sec%transport( 3,jclass) = 0._wp
945	sec%transport( 4,jclass) = 0._wp
946	sec%transport( 5,jclass) = 0._wp
947	sec%transport( 6,jclass) = 0._wp
948	ENDIF
949
950	ENDIF ! end of test if point is in class
951
952	END DO ! end of loop on the classes
953
954	END DO ! loop over jk
955
956	#if defined key_si3
957
958	!ICE CASE
959	IF( sec%ll_ice_section )THEN
960
961	IF ( transports_2d(1,jsec,jseg) .GE. 0.0 ) THEN
962	sec%transport( 7,1) = sec%transport( 7,1)+transports_2d(1,jsec,jseg)*1.E-6
963	ELSE
964	sec%transport( 8,1) = sec%transport( 8,1)+transports_2d(1,jsec,jseg)*1.E-6
965	ENDIF
966
967	IF ( transports_2d(3,jsec,jseg) .GE. 0.0 ) THEN
968	sec%transport( 9,1) = sec%transport( 9,1)+transports_2d(2,jsec,jseg)*1.E-6
969	ELSE
970	sec%transport(10,1) = sec%transport(10,1)+transports_2d(2,jsec,jseg)*1.E-6
971	ENDIF
972
973	ENDIF !end of ice case
974	#endif
975
976	END DO !end of loop on the segment
977
978	ELSE !if sec%nb_point =0
979	sec%transport(1:2,:)=0.
980	IF (sec%llstrpond) sec%transport(3:6,:)=0.
981	IF (sec%ll_ice_section) sec%transport(7:10,:)=0.
982	ENDIF !end of sec%nb_point =0 case
983
984	END SUBROUTINE dia_dct_sum
985
986
987	SUBROUTINE dia_dct_wri(kt,ksec,sec)
988	!!-------------------------------------------------------------
989	!! Write transport output in numdct
990	!!
991	!! Purpose: Write transports in ascii files
992	!!
993	!! Method:
994	!! 1. Write volume transports in "volume_transport"
995	!! Unit: Sv : area * Velocity / 1.e6
996	!!
997	!! 2. Write heat transports in "heat_transport"
998	!! Unit: Peta W : area * Velocity * T * rhop * Cp * 1.e-15
999	!!
1000	!! 3. Write salt transports in "salt_transport"
1001	!! Unit: 10^9 Kg/m^2/s : area * Velocity * S * rhop * 1.e-9
1002	!!
1003	!!-------------------------------------------------------------
1004	!!arguments
1005	INTEGER, INTENT(IN) :: kt ! time-step
1006	TYPE(SECTION), INTENT(INOUT) :: sec ! section to write
1007	INTEGER ,INTENT(IN) :: ksec ! section number
1008
1009	!!local declarations
1010	INTEGER :: jclass ! Dummy loop
1011	CHARACTER(len=2) :: classe ! Classname
1012	REAL(wp) :: zbnd1,zbnd2 ! Class bounds
1013	REAL(wp) :: zslope ! section's slope coeff
1014	!
1015	REAL(wp), DIMENSION(nb_type_class):: zsumclasses ! 1D workspace
1016	!!-------------------------------------------------------------
1017
1018	zsumclasses(:)=0._wp
1019	zslope = sec%slopeSection
1020
1021
1022	DO jclass=1,MAX(1,sec%nb_class-1)
1023
1024	classe = 'N '
1025	zbnd1 = 0._wp
1026	zbnd2 = 0._wp
1027	zsumclasses(1:nb_type_class)=zsumclasses(1:nb_type_class)+sec%transport(1:nb_type_class,jclass)
1028
1029
1030	!insitu density classes transports
1031	IF( ( sec%zsigi(jclass) .NE. 99._wp ) .AND. &
1032	( sec%zsigi(jclass+1) .NE. 99._wp ) )THEN
1033	classe = 'DI '
1034	zbnd1 = sec%zsigi(jclass)
1035	zbnd2 = sec%zsigi(jclass+1)
1036	ENDIF
1037	!potential density classes transports
1038	IF( ( sec%zsigp(jclass) .NE. 99._wp ) .AND. &
1039	( sec%zsigp(jclass+1) .NE. 99._wp ) )THEN
1040	classe = 'DP '
1041	zbnd1 = sec%zsigp(jclass)
1042	zbnd2 = sec%zsigp(jclass+1)
1043	ENDIF
1044	!depth classes transports
1045	IF( ( sec%zlay(jclass) .NE. 99._wp ) .AND. &
1046	( sec%zlay(jclass+1) .NE. 99._wp ) )THEN
1047	classe = 'Z '
1048	zbnd1 = sec%zlay(jclass)
1049	zbnd2 = sec%zlay(jclass+1)
1050	ENDIF
1051	!salinity classes transports
1052	IF( ( sec%zsal(jclass) .NE. 99._wp ) .AND. &
1053	( sec%zsal(jclass+1) .NE. 99._wp ) )THEN
1054	classe = 'S '
1055	zbnd1 = sec%zsal(jclass)
1056	zbnd2 = sec%zsal(jclass+1)
1057	ENDIF
1058	!temperature classes transports
1059	IF( ( sec%ztem(jclass) .NE. 99._wp ) .AND. &
1060	( sec%ztem(jclass+1) .NE. 99._wp ) ) THEN
1061	classe = 'T '
1062	zbnd1 = sec%ztem(jclass)
1063	zbnd2 = sec%ztem(jclass+1)
1064	ENDIF
1065
1066	!write volume transport per class
1067	WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, &
1068	jclass,classe,zbnd1,zbnd2,&
1069	sec%transport(1,jclass),sec%transport(2,jclass), &
1070	sec%transport(1,jclass)+sec%transport(2,jclass)
1071
1072	IF( sec%llstrpond )THEN
1073
1074	!write heat transport per class:
1075	WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope, &
1076	jclass,classe,zbnd1,zbnd2,&
1077	sec%transport(3,jclass)1.e-15,sec%transport(4,jclass)1.e-15, &
1078	( sec%transport(3,jclass)+sec%transport(4,jclass) )*1.e-15
1079	!write salt transport per class
1080	WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope, &
1081	jclass,classe,zbnd1,zbnd2,&
1082	sec%transport(5,jclass)1.e-9,sec%transport(6,jclass)1.e-9,&
1083	(sec%transport(5,jclass)+sec%transport(6,jclass))*1.e-9
1084	ENDIF
1085
1086	ENDDO
1087
1088	zbnd1 = 0._wp
1089	zbnd2 = 0._wp
1090	jclass=0
1091
1092	!write total volume transport
1093	WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, &
1094	jclass,"total",zbnd1,zbnd2,&
1095	zsumclasses(1),zsumclasses(2),zsumclasses(1)+zsumclasses(2)
1096
1097	IF( sec%llstrpond )THEN
1098
1099	!write total heat transport
1100	WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope, &
1101	jclass,"total",zbnd1,zbnd2,&
1102	zsumclasses(3)1.e-15,zsumclasses(4)1.e-15,&
1103	(zsumclasses(3)+zsumclasses(4) )*1.e-15
1104	!write total salt transport
1105	WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope, &
1106	jclass,"total",zbnd1,zbnd2,&
1107	zsumclasses(5)1.e-9,zsumclasses(6)1.e-9,&
1108	(zsumclasses(5)+zsumclasses(6))*1.e-9
1109	ENDIF
1110
1111
1112	IF ( sec%ll_ice_section) THEN
1113	!write total ice volume transport
1114	WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,&
1115	jclass,"ice_vol",zbnd1,zbnd2,&
1116	sec%transport(7,1),sec%transport(8,1),&
1117	sec%transport(7,1)+sec%transport(8,1)
1118	!write total ice surface transport
1119	WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,&
1120	jclass,"ice_surf",zbnd1,zbnd2,&
1121	sec%transport(9,1),sec%transport(10,1), &
1122	sec%transport(9,1)+sec%transport(10,1)
1123	ENDIF
1124
1125	118 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3F12.4)
1126	119 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6)
1127	!
1128	END SUBROUTINE dia_dct_wri
1129
1130
1131	FUNCTION interp(ki, kj, kk, cd_point, ptab)
1132	!!----------------------------------------------------------------------
1133	!!
1134	!! Purpose: compute temperature/salinity/density at U-point or V-point
1135	!! --------
1136	!!
1137	!! Method:
1138	!! ------
1139	!!
1140	!! ====> full step and partial step
1141	!!
1142	!!
1143	!! \| I \| I+1 \| Z=temperature/salinity/density at U-poinT
1144	!! \| \| \|
1145	!! ---------------------------------------- 1. Veritcal interpolation: compute zbis
1146	!! \| \| \| interpolation between ptab(I,J,K) and ptab(I,J,K+1)
1147	!! \| \| \| zbis =
1148	!! \| \| \| [ e3w(I+1,J,K)ptab(I,J,K) + ( e3w(I,J,K) - e3w(I+1,J,K) ) ptab(I,J,K-1) ]
1149	!! \| \| \| /[ e3w(I+1,J,K) + e3w(I,J,K) - e3w(I+1,J,K) ]
1150	!! \| \| \|
1151	!! \| \| \| 2. Horizontal interpolation: compute value at U/V point
1152	!!K-1 \| ptab(I,J,K-1) \| \| interpolation between zbis and ptab(I+1,J,K)
1153	!! \| . \| \|
1154	!! \| . \| \| interp = ( 0.5zet2zbis + 0.5zet1ptab(I+1,J,K) )/(0.5zet2+0.5zet1)
1155	!! \| . \| \|
1156	!! ------------------------------------------
1157	!! \| . \| \|
1158	!! \| . \| \|
1159	!! \| . \| \|
1160	!!K \| zbis.......U...ptab(I+1,J,K) \|
1161	!! \| . \| \|
1162	!! \| ptab(I,J,K) \| \|
1163	!! \| \|------------------\|
1164	!! \| \| partials \|
1165	!! \| \| steps \|
1166	!! -------------------------------------------
1167	!! <----zet1------><----zet2--------->
1168	!!
1169	!!
1170	!! ====> s-coordinate
1171	!!
1172	!! \| \| \| 1. Compute distance between T1 and U points: SQRT( zdep1^2 + (0.5 * zet1 )^2
1173	!! \| \| \| Compute distance between T2 and U points: SQRT( zdep2^2 + (0.5 * zet2 )^2
1174	!! \| \| ptab(I+1,J,K) \|
1175	!! \| \| T2 \| 2. Interpolation between T1 and T2 values at U point
1176	!! \| \| ^ \|
1177	!! \| \| \| zdep2 \|
1178	!! \| \| \| \|
1179	!! \| ^ U v \|
1180	!! \| \| \| \|
1181	!! \| \| zdep1 \| \|
1182	!! \| v \| \|
1183	!! \| T1 \| \|
1184	!! \| ptab(I,J,K) \| \|
1185	!! \| \| \|
1186	!! \| \| \|
1187	!!
1188	!! <----zet1--------><----zet2--------->
1189	!!
1190	!!----------------------------------------------------------------------
1191	!*arguments
1192	INTEGER, INTENT(IN) :: ki, kj, kk ! coordinate of point
1193	CHARACTER(len=1), INTENT(IN) :: cd_point ! type of point (U, V)
1194	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ptab ! variable to compute at (ki, kj, kk )
1195	REAL(wp) :: interp ! interpolated variable
1196
1197	!*local declations
1198	INTEGER :: ii1, ij1, ii2, ij2 ! local integer
1199	REAL(wp):: ze3t, ze3, zwgt1, zwgt2, zbis, zdepu ! local real
1200	REAL(wp):: zet1, zet2 ! weight for interpolation
1201	REAL(wp):: zdep1,zdep2 ! differences of depth
1202	REAL(wp):: zmsk ! mask value
1203	!!----------------------------------------------------------------------
1204
1205	IF( cd_point=='U' )THEN
1206	ii1 = ki ; ij1 = kj
1207	ii2 = ki+1 ; ij2 = kj
1208
1209	zet1=e1t(ii1,ij1)
1210	zet2=e1t(ii2,ij2)
1211	zmsk=umask(ii1,ij1,kk)
1212
1213
1214	ELSE ! cd_point=='V'
1215	ii1 = ki ; ij1 = kj
1216	ii2 = ki ; ij2 = kj+1
1217
1218	zet1=e2t(ii1,ij1)
1219	zet2=e2t(ii2,ij2)
1220	zmsk=vmask(ii1,ij1,kk)
1221
1222	ENDIF
1223
1224	IF( ln_sco )THEN ! s-coordinate case
1225
1226	zdepu = ( gdept_n(ii1,ij1,kk) + gdept_n(ii2,ij2,kk) ) * 0.5_wp
1227	zdep1 = gdept_n(ii1,ij1,kk) - zdepu
1228	zdep2 = gdept_n(ii2,ij2,kk) - zdepu
1229
1230	! weights
1231	zwgt1 = SQRT( ( 0.5 * zet1 ) * ( 0.5 * zet1 ) + ( zdep1 * zdep1 ) )
1232	zwgt2 = SQRT( ( 0.5 * zet2 ) * ( 0.5 * zet2 ) + ( zdep2 * zdep2 ) )
1233
1234	! result
1235	interp = zmsk * ( zwgt2 * ptab(ii1,ij1,kk) + zwgt1 * ptab(ii1,ij1,kk) ) / ( zwgt2 + zwgt1 )
1236
1237
1238	ELSE ! full step or partial step case
1239
1240	ze3t = e3t_n(ii2,ij2,kk) - e3t_n(ii1,ij1,kk)
1241	zwgt1 = ( e3w_n(ii2,ij2,kk) - e3w_n(ii1,ij1,kk) ) / e3w_n(ii2,ij2,kk)
1242	zwgt2 = ( e3w_n(ii1,ij1,kk) - e3w_n(ii2,ij2,kk) ) / e3w_n(ii1,ij1,kk)
1243
1244	IF(kk .NE. 1)THEN
1245
1246	IF( ze3t >= 0. )THEN
1247	! zbis
1248	zbis = ptab(ii2,ij2,kk) + zwgt1 * ( ptab(ii2,ij2,kk-1) - ptab(ii2,ij2,kk) )
1249	! result
1250	interp = zmsk * ( zet2 * ptab(ii1,ij1,kk) + zet1 * zbis )/( zet1 + zet2 )
1251	ELSE
1252	! zbis
1253	zbis = ptab(ii1,ij1,kk) + zwgt2 * ( ptab(ii1,ij1,kk-1) - ptab(ii1,ij2,kk) )
1254	! result
1255	interp = zmsk * ( zet2 * zbis + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 )
1256	ENDIF
1257
1258	ELSE
1259	interp = zmsk * ( zet2 * ptab(ii1,ij1,kk) + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 )
1260	ENDIF
1261
1262	ENDIF
1263	!
1264	END FUNCTION interp
1265
1266	#else
1267	!!----------------------------------------------------------------------
1268	!! Default option : Dummy module
1269	!!----------------------------------------------------------------------
1270	LOGICAL, PUBLIC, PARAMETER :: lk_diadct = .FALSE. !: diamht flag
1271	PUBLIC
1272	!! $Id$
1273	CONTAINS
1274
1275	SUBROUTINE dia_dct_init ! Dummy routine
1276	IMPLICIT NONE
1277	WRITE(,) 'dia_dct_init: You should not have seen this print! error?'
1278	END SUBROUTINE dia_dct_init
1279
1280	SUBROUTINE dia_dct( kt ) ! Dummy routine
1281	IMPLICIT NONE
1282	INTEGER, INTENT( in ) :: kt ! ocean time-step index
1283	WRITE(,) 'dia_dct: You should not have seen this print! error?', kt
1284	END SUBROUTINE dia_dct
1285	#endif
1286
1287	!!======================================================================
1288	END MODULE diadct

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

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

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0_mirror_text_diagnostics/src/OCE/DIA/diadct.F90 @ 10986

Download in other formats: