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/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: branches/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 9366

Last change on this file since 9366 was 9366, checked in by andmirek, 6 years ago
#2050 first version. Compiled OK in moci test suite
Property svn:executable set to ``*
File size: 60.5 KB

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

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

Download in other formats: