Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

diadct.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago
All wrk_alloc removed
Property svn:keywords set to `Id`
File size: 59.3 KB

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

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

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

Context Navigation

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 7910

Download in other formats: