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

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

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

Context Navigation

source: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 8058

Last change on this file since 8058 was 8058, checked in by jgraham, 7 years ago
Clear keywords
Property svn:executable set to ``*
File size: 59.9 KB

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

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

Download in other formats: