Context Navigation

diadct.F90 @ 8563

Last change on this file since 8563 was 8563, checked in by clem, 7 years ago
change variable names (ht_s => h_s & ht_i => h_i)
Property svn:keywords set to `Id`
File size: 59.2 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_lim3
35	USE ice
36	#endif
37	USE domvvl
38	USE timing ! preformance summary
39	USE wrk_nemo ! working arrays
40
41	IMPLICIT NONE
42	PRIVATE
43
44	PUBLIC dia_dct ! routine called by step.F90
45	PUBLIC dia_dct_init ! routine called by opa.F90
46	PUBLIC diadct_alloc ! routine called by nemo_init in nemogcm.F90
47	PRIVATE readsec
48	PRIVATE removepoints
49	PRIVATE transport
50	PRIVATE dia_dct_wri
51
52	LOGICAL, PUBLIC, PARAMETER :: lk_diadct = .TRUE. !: model-data diagnostics flag
53
54	INTEGER :: nn_dct ! Frequency of computation
55	INTEGER :: nn_dctwri ! Frequency of output
56	INTEGER :: nn_secdebug ! Number of the section to debug
57
58	INTEGER, PARAMETER :: nb_class_max = 10
59	INTEGER, PARAMETER :: nb_sec_max = 150
60	INTEGER, PARAMETER :: nb_point_max = 2000
61	INTEGER, PARAMETER :: nb_type_class = 10
62	INTEGER, PARAMETER :: nb_3d_vars = 3
63	INTEGER, PARAMETER :: nb_2d_vars = 2
64	INTEGER :: nb_sec
65
66	TYPE POINT_SECTION
67	INTEGER :: I,J
68	END TYPE POINT_SECTION
69
70	TYPE COORD_SECTION
71	REAL(wp) :: lon,lat
72	END TYPE COORD_SECTION
73
74	TYPE SECTION
75	CHARACTER(len=60) :: name ! name of the sec
76	LOGICAL :: llstrpond ! true if you want the computation of salt and
77	! heat transports
78	LOGICAL :: ll_ice_section ! ice surface and ice volume computation
79	LOGICAL :: ll_date_line ! = T if the section crosses the date-line
80	TYPE(COORD_SECTION), DIMENSION(2) :: coordSec ! longitude and latitude of the extremities of the sec
81	INTEGER :: nb_class ! number of boundaries for density classes
82	INTEGER, DIMENSION(nb_point_max) :: direction ! vector direction of the point in the section
83	CHARACTER(len=40),DIMENSION(nb_class_max) :: classname ! characteristics of the class
84	REAL(wp), DIMENSION(nb_class_max) :: zsigi ,&! in-situ density classes (99 if you don't want)
85	zsigp ,&! potential density classes (99 if you don't want)
86	zsal ,&! salinity classes (99 if you don't want)
87	ztem ,&! temperature classes(99 if you don't want)
88	zlay ! level classes (99 if you don't want)
89	REAL(wp), DIMENSION(nb_type_class,nb_class_max) :: transport ! transport output
90	REAL(wp) :: slopeSection ! slope of the section
91	INTEGER :: nb_point ! number of points in the section
92	TYPE(POINT_SECTION),DIMENSION(nb_point_max) :: listPoint ! list of points in the sections
93	END TYPE SECTION
94
95	TYPE(SECTION),DIMENSION(nb_sec_max) :: secs ! Array of sections
96
97	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: transports_3d
98	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: transports_2d
99
100	!!----------------------------------------------------------------------
101	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
102	!! $Id$
103	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
104	!!----------------------------------------------------------------------
105	CONTAINS
106
107	INTEGER FUNCTION diadct_alloc()
108	!!----------------------------------------------------------------------
109	!! * FUNCTION diadct_alloc *
110	!!----------------------------------------------------------------------
111	INTEGER :: ierr(2)
112	!!----------------------------------------------------------------------
113
114	ALLOCATE(transports_3d(nb_3d_vars,nb_sec_max,nb_point_max,jpk), STAT=ierr(1) )
115	ALLOCATE(transports_2d(nb_2d_vars,nb_sec_max,nb_point_max) , STAT=ierr(2) )
116
117	diadct_alloc = MAXVAL( ierr )
118	IF( diadct_alloc /= 0 ) CALL ctl_warn('diadct_alloc: failed to allocate arrays')
119
120	END FUNCTION diadct_alloc
121
122
123	SUBROUTINE dia_dct_init
124	!!---------------------------------------------------------------------
125	!! * ROUTINE diadct *
126	!!
127	!! ** Purpose: Read the namelist parameters
128	!! Open output files
129	!!
130	!!---------------------------------------------------------------------
131	NAMELIST/namdct/nn_dct,nn_dctwri,nn_secdebug
132	INTEGER :: ios ! Local integer output status for namelist read
133
134	IF( nn_timing == 1 ) CALL timing_start('dia_dct_init')
135
136	REWIND( numnam_ref ) ! Namelist namdct in reference namelist : Diagnostic: transport through sections
137	READ ( numnam_ref, namdct, IOSTAT = ios, ERR = 901)
138	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdct in reference namelist', lwp )
139
140	REWIND( numnam_cfg ) ! Namelist namdct in configuration namelist : Diagnostic: transport through sections
141	READ ( numnam_cfg, namdct, IOSTAT = ios, ERR = 902 )
142	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdct in configuration namelist', lwp )
143	IF(lwm) WRITE ( numond, namdct )
144
145	IF( lwp ) THEN
146	WRITE(numout,*) " "
147	WRITE(numout,*) "diadct_init: compute transports through sections "
148	WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
149	WRITE(numout,*) " Frequency of computation: nn_dct = ",nn_dct
150	WRITE(numout,*) " Frequency of write: nn_dctwri = ",nn_dctwri
151
152	IF ( nn_secdebug .GE. 1 .AND. nn_secdebug .LE. nb_sec_max )THEN
153	WRITE(numout,*)" Debug section number: ", nn_secdebug
154	ELSE IF ( nn_secdebug == 0 )THEN ; WRITE(numout,*)" No section to debug"
155	ELSE IF ( nn_secdebug == -1 )THEN ; WRITE(numout,*)" Debug all sections"
156	ELSE ; WRITE(numout,*)" Wrong value for nn_secdebug : ",nn_secdebug
157	ENDIF
158
159	IF(nn_dct .GE. nn_dctwri .AND. MOD(nn_dct,nn_dctwri) .NE. 0) &
160	& CALL ctl_stop( 'diadct: nn_dct should be smaller and a multiple of nn_dctwri' )
161
162	ENDIF
163
164	!Read section_ijglobal.diadct
165	CALL readsec
166
167	!open output file
168	IF( lwm ) THEN
169	CALL ctl_opn( numdct_vol, 'volume_transport', 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
170	CALL ctl_opn( numdct_heat, 'heat_transport' , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
171	CALL ctl_opn( numdct_salt, 'salt_transport' , 'NEW', 'FORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
172	ENDIF
173
174	! Initialise arrays to zero
175	transports_3d(:,:,:,:)=0.0
176	transports_2d(:,:,:) =0.0
177
178	IF( nn_timing == 1 ) CALL timing_stop('dia_dct_init')
179	!
180	END SUBROUTINE dia_dct_init
181
182
183	SUBROUTINE dia_dct(kt)
184	!!---------------------------------------------------------------------
185	!! * ROUTINE diadct *
186	!!
187	!! Purpose :: Compute section transports and write it in numdct files
188	!!
189	!! Method :: All arrays initialised to zero in dct_init
190	!! Each nn_dct time step call subroutine 'transports' for
191	!! each section to sum the transports over each grid cell.
192	!! Each nn_dctwri time step:
193	!! Divide the arrays by the number of summations to gain
194	!! an average value
195	!! Call dia_dct_sum to sum relevant grid boxes to obtain
196	!! totals for each class (density, depth, temp or sal)
197	!! Call dia_dct_wri to write the transports into file
198	!! Reinitialise all relevant arrays to zero
199	!!---------------------------------------------------------------------
200	INTEGER,INTENT(in) ::kt
201	!
202	INTEGER :: jsec, &! loop on sections
203	itotal ! nb_sec_maxnb_type_classnb_class_max
204	LOGICAL :: lldebug =.FALSE. ! debug a section
205
206	INTEGER , DIMENSION(1) :: ish ! tmp array for mpp_sum
207	INTEGER , DIMENSION(3) :: ish2 ! "
208	REAL(wp), POINTER, DIMENSION(:) :: zwork ! "
209	REAL(wp), POINTER, DIMENSION(:,:,:):: zsum ! "
210	!!---------------------------------------------------------------------
211	!
212	IF( nn_timing == 1 ) CALL timing_start('dia_dct')
213
214	IF( lk_mpp )THEN
215	itotal = nb_sec_maxnb_type_classnb_class_max
216	CALL wrk_alloc( itotal , zwork )
217	CALL wrk_alloc( nb_sec_max,nb_type_class,nb_class_max , zsum )
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 ) 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( 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	CALL wrk_dealloc( itotal , zwork )
289	CALL wrk_dealloc( nb_sec_max,nb_type_class,nb_class_max , zsum )
290	ENDIF
291
292	IF( nn_timing == 1 ) CALL timing_stop('dia_dct')
293	!
294	END SUBROUTINE dia_dct
295
296	SUBROUTINE readsec
297	!!---------------------------------------------------------------------
298	!! * ROUTINE readsec *
299	!!
300	!! ** Purpose:
301	!! Read a binary file(section_ijglobal.diadct)
302	!! generated by the tools "NEMOGCM/TOOLS/SECTIONS_DIADCT"
303	!!
304	!!
305	!!---------------------------------------------------------------------
306	!! * Local variables
307	INTEGER :: iptglo , iptloc ! Global and local number of points for a section
308	INTEGER :: isec, iiglo, ijglo, iiloc, ijloc,iost,i1 ,i2 ! temporary integer
309	INTEGER :: jsec, jpt ! dummy loop indices
310
311	INTEGER, DIMENSION(2) :: icoord
312	CHARACTER(len=160) :: clname !filename
313	CHARACTER(len=200) :: cltmp
314	CHARACTER(len=200) :: clformat !automatic format
315	TYPE(POINT_SECTION),DIMENSION(nb_point_max) ::coordtemp !contains listpoints coordinates
316	!read in the file
317	INTEGER, POINTER, DIMENSION(:) :: directemp !contains listpoints directions
318	!read in the files
319	LOGICAL :: llbon ,&!local logical
320	lldebug !debug the section
321	!!-------------------------------------------------------------------------------------
322	CALL wrk_alloc( nb_point_max, directemp )
323
324	!open input file
325	!---------------
326	CALL ctl_opn( numdct_in, 'section_ijglobal.diadct', 'OLD', 'UNFORMATTED', 'SEQUENTIAL', -1, numout, .FALSE. )
327
328	!---------------
329	!Read input file
330	!---------------
331
332	DO jsec=1,nb_sec_max !loop on the nb_sec sections
333
334	IF ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) &
335	& WRITE(numout,*)'debuging for section number: ',jsec
336
337	!initialization
338	!---------------
339	secs(jsec)%name=''
340	secs(jsec)%llstrpond = .FALSE. ; secs(jsec)%ll_ice_section = .FALSE.
341	secs(jsec)%ll_date_line = .FALSE. ; secs(jsec)%nb_class = 0
342	secs(jsec)%zsigi = 99._wp ; secs(jsec)%zsigp = 99._wp
343	secs(jsec)%zsal = 99._wp ; secs(jsec)%ztem = 99._wp
344	secs(jsec)%zlay = 99._wp
345	secs(jsec)%transport = 0._wp ; secs(jsec)%nb_point = 0
346
347	!read section's number / name / computing choices / classes / slopeSection / points number
348	!-----------------------------------------------------------------------------------------
349	READ(numdct_in,iostat=iost)isec
350	IF (iost .NE. 0 )EXIT !end of file
351	WRITE(cltmp,'(a,i4.4,a,i4.4)')'diadct: read sections : Problem of section number: isec= ',isec,' and jsec= ',jsec
352	IF( jsec .NE. isec ) CALL ctl_stop( cltmp )
353
354	IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )WRITE(numout,*)"isec ",isec
355
356	READ(numdct_in)secs(jsec)%name
357	READ(numdct_in)secs(jsec)%llstrpond
358	READ(numdct_in)secs(jsec)%ll_ice_section
359	READ(numdct_in)secs(jsec)%ll_date_line
360	READ(numdct_in)secs(jsec)%coordSec
361	READ(numdct_in)secs(jsec)%nb_class
362	READ(numdct_in)secs(jsec)%zsigi
363	READ(numdct_in)secs(jsec)%zsigp
364	READ(numdct_in)secs(jsec)%zsal
365	READ(numdct_in)secs(jsec)%ztem
366	READ(numdct_in)secs(jsec)%zlay
367	READ(numdct_in)secs(jsec)%slopeSection
368	READ(numdct_in)iptglo
369
370	!debug
371	!-----
372
373	IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )THEN
374
375	WRITE(clformat,'(a,i2,a)') '(A40,', nb_class_max,'(f8.3,1X))'
376
377	WRITE(numout,*) " Section name : ",TRIM(secs(jsec)%name)
378	WRITE(numout,*) " Compute heat and salt transport ? ",secs(jsec)%llstrpond
379	WRITE(numout,*) " Compute ice transport ? ",secs(jsec)%ll_ice_section
380	WRITE(numout,*) " Section crosses date-line ? ",secs(jsec)%ll_date_line
381	WRITE(numout,*) " Slope section : ",secs(jsec)%slopeSection
382	WRITE(numout,*) " Number of points in the section: ",iptglo
383	WRITE(numout,*) " Number of classes ",secs(jsec)%nb_class
384	WRITE(numout,clformat)" Insitu density classes : ",secs(jsec)%zsigi
385	WRITE(numout,clformat)" Potential density classes : ",secs(jsec)%zsigp
386	WRITE(numout,clformat)" Salinity classes : ",secs(jsec)%zsal
387	WRITE(numout,clformat)" Temperature classes : ",secs(jsec)%ztem
388	WRITE(numout,clformat)" Depth classes : ",secs(jsec)%zlay
389	ENDIF
390
391	IF( iptglo /= 0 )THEN
392
393	!read points'coordinates and directions
394	!--------------------------------------
395	coordtemp(:) = POINT_SECTION(0,0) !list of points read
396	directemp(:) = 0 !value of directions of each points
397	DO jpt=1,iptglo
398	READ(numdct_in) i1, i2
399	coordtemp(jpt)%I = i1
400	coordtemp(jpt)%J = i2
401	ENDDO
402	READ(numdct_in) directemp(1:iptglo)
403
404	!debug
405	!-----
406	IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )THEN
407	WRITE(numout,*)" List of points in global domain:"
408	DO jpt=1,iptglo
409	WRITE(numout,*)' # I J ',jpt,coordtemp(jpt),directemp(jpt)
410	ENDDO
411	ENDIF
412
413	!Now each proc selects only points that are in its domain:
414	!--------------------------------------------------------
415	iptloc = 0 ! initialize number of points selected
416	DO jpt = 1, iptglo ! loop on listpoint read in the file
417	!
418	iiglo=coordtemp(jpt)%I ! global coordinates of the point
419	ijglo=coordtemp(jpt)%J ! "
420
421	IF( iiglo==jpiglo .AND. nimpp==1 ) iiglo = 2 !!gm BUG: Hard coded periodicity !
422
423	iiloc=iiglo-nimpp+1 ! local coordinates of the point
424	ijloc=ijglo-njmpp+1 ! "
425
426	!verify if the point is on the local domain:(1,nlei)*(1,nlej)
427	IF( iiloc >= 1 .AND. iiloc <= nlei .AND. &
428	ijloc >= 1 .AND. ijloc <= nlej )THEN
429	iptloc = iptloc + 1 ! count local points
430	secs(jsec)%listPoint(iptloc) = POINT_SECTION(mi0(iiglo),mj0(ijglo)) ! store local coordinates
431	secs(jsec)%direction(iptloc) = directemp(jpt) ! store local direction
432	ENDIF
433	!
434	END DO
435
436	secs(jsec)%nb_point=iptloc !store number of section's points
437
438	!debug
439	!-----
440	IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )THEN
441	WRITE(numout,*)" List of points selected by the proc:"
442	DO jpt = 1,iptloc
443	iiglo = secs(jsec)%listPoint(jpt)%I + nimpp - 1
444	ijglo = secs(jsec)%listPoint(jpt)%J + njmpp - 1
445	WRITE(numout,*)' # I J : ',iiglo,ijglo
446	ENDDO
447	ENDIF
448
449	IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN
450	DO jpt = 1,iptloc
451	iiglo = secs(jsec)%listPoint(jpt)%I + nimpp - 1
452	ijglo = secs(jsec)%listPoint(jpt)%J + njmpp - 1
453	ENDDO
454	ENDIF
455
456	!remove redundant points between processors
457	!------------------------------------------
458	lldebug = .FALSE. ; IF ( jsec==nn_secdebug .OR. nn_secdebug==-1 ) lldebug = .TRUE.
459	IF( iptloc .NE. 0 )THEN
460	CALL removepoints(secs(jsec),'I','top_list',lldebug)
461	CALL removepoints(secs(jsec),'I','bot_list',lldebug)
462	CALL removepoints(secs(jsec),'J','top_list',lldebug)
463	CALL removepoints(secs(jsec),'J','bot_list',lldebug)
464	ENDIF
465	IF(jsec==nn_secdebug .AND. secs(jsec)%nb_point .NE. 0)THEN
466	DO jpt = 1,secs(jsec)%nb_point
467	iiglo = secs(jsec)%listPoint(jpt)%I + nimpp - 1
468	ijglo = secs(jsec)%listPoint(jpt)%J + njmpp - 1
469	ENDDO
470	ENDIF
471
472	!debug
473	!-----
474	IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )THEN
475	WRITE(numout,*)" List of points after removepoints:"
476	iptloc = secs(jsec)%nb_point
477	DO jpt = 1,iptloc
478	iiglo = secs(jsec)%listPoint(jpt)%I + nimpp - 1
479	ijglo = secs(jsec)%listPoint(jpt)%J + njmpp - 1
480	WRITE(numout,*)' # I J : ',iiglo,ijglo
481	CALL FLUSH(numout)
482	ENDDO
483	ENDIF
484
485	ELSE ! iptglo = 0
486	IF( jsec==nn_secdebug .OR. nn_secdebug==-1 )&
487	WRITE(numout,*)' No points for this section.'
488	ENDIF
489
490	ENDDO !end of the loop on jsec
491
492	nb_sec = jsec-1 !number of section read in the file
493
494	CALL wrk_dealloc( nb_point_max, directemp )
495	!
496	END SUBROUTINE readsec
497
498	SUBROUTINE removepoints(sec,cdind,cdextr,ld_debug)
499	!!---------------------------------------------------------------------------
500	!! *** function removepoints
501	!!
502	!! ** Purpose :: Remove points which are common to 2 procs
503	!!
504	!----------------------------------------------------------------------------
505	!! * arguments
506	TYPE(SECTION),INTENT(INOUT) :: sec
507	CHARACTER(len=1),INTENT(IN) :: cdind ! = 'I'/'J'
508	CHARACTER(len=8),INTENT(IN) :: cdextr ! = 'top_list'/'bot_list'
509	LOGICAL,INTENT(IN) :: ld_debug
510
511	!! * Local variables
512	INTEGER :: iextr ,& !extremity of listpoint that we verify
513	iind ,& !coord of listpoint that we verify
514	itest ,& !indice value of the side of the domain
515	!where points could be redundant
516	isgn ,& ! isgn= 1 : scan listpoint from start to end
517	! isgn=-1 : scan listpoint from end to start
518	istart,iend !first and last points selected in listpoint
519	INTEGER :: jpoint !loop on list points
520	INTEGER, POINTER, DIMENSION(:) :: idirec !contains temporary sec%direction
521	INTEGER, POINTER, DIMENSION(:,:) :: icoord !contains temporary sec%listpoint
522	!----------------------------------------------------------------------------
523	CALL wrk_alloc( nb_point_max, idirec )
524	CALL wrk_alloc( 2, nb_point_max, icoord )
525
526	IF( ld_debug )WRITE(numout,*)' -------------------------'
527	IF( ld_debug )WRITE(numout,*)' removepoints in listpoint'
528
529	!iextr=extremity of list_point that we verify
530	IF ( cdextr=='bot_list' )THEN ; iextr=1 ; isgn=1
531	ELSE IF ( cdextr=='top_list' )THEN ; iextr=sec%nb_point ; isgn=-1
532	ELSE ; CALL ctl_stop("removepoints :Wrong value for cdextr")
533	ENDIF
534
535	!which coordinate shall we verify ?
536	IF ( cdind=='I' )THEN ; itest=nlei ; iind=1
537	ELSE IF ( cdind=='J' )THEN ; itest=nlej ; iind=2
538	ELSE ; CALL ctl_stop("removepoints :Wrong value for cdind")
539	ENDIF
540
541	IF( ld_debug )THEN
542	WRITE(numout,*)' case: coord/list extr/domain side'
543	WRITE(numout,*)' ', cdind,' ',cdextr,' ',itest
544	WRITE(numout,*)' Actual number of points: ',sec%nb_point
545	ENDIF
546
547	icoord(1,1:nb_point_max) = sec%listPoint%I
548	icoord(2,1:nb_point_max) = sec%listPoint%J
549	idirec = sec%direction
550	sec%listPoint = POINT_SECTION(0,0)
551	sec%direction = 0
552
553	jpoint=iextr+isgn
554	DO WHILE( jpoint .GE. 1 .AND. jpoint .LE. sec%nb_point )
555	IF( icoord( iind,jpoint-isgn ) == itest .AND. icoord( iind,jpoint ) == itest )THEN ; jpoint=jpoint+isgn
556	ELSE ; EXIT
557	ENDIF
558	ENDDO
559
560	IF( cdextr=='bot_list')THEN ; istart=jpoint-1 ; iend=sec%nb_point
561	ELSE ; istart=1 ; iend=jpoint+1
562	ENDIF
563
564	sec%listPoint(1:1+iend-istart)%I = icoord(1,istart:iend)
565	sec%listPoint(1:1+iend-istart)%J = icoord(2,istart:iend)
566	sec%direction(1:1+iend-istart) = idirec(istart:iend)
567	sec%nb_point = iend-istart+1
568
569	IF( ld_debug )THEN
570	WRITE(numout,*)' Number of points after removepoints :',sec%nb_point
571	WRITE(numout,*)' sec%direction after removepoints :',sec%direction(1:sec%nb_point)
572	ENDIF
573
574	CALL wrk_dealloc( nb_point_max, idirec )
575	CALL wrk_dealloc( 2, nb_point_max, icoord )
576	END SUBROUTINE removepoints
577
578	SUBROUTINE transport(sec,ld_debug,jsec)
579	!!-------------------------------------------------------------------------------------------
580	!! * ROUTINE transport *
581	!!
582	!! Purpose :: Compute the transport for each point in a section
583	!!
584	!! Method :: Loop over each segment, and each vertical level and add the transport
585	!! Be aware :
586	!! One section is a sum of segments
587	!! One segment is defined by 2 consecutive points in sec%listPoint
588	!! All points of sec%listPoint are positioned on the F-point of the cell
589	!!
590	!! There are two loops:
591	!! loop on the segment between 2 nodes
592	!! loop on the level jk !!
593	!!
594	!! Output :: Arrays containing the volume,density,heat,salt transports for each i
595	!! point in a section, summed over each nn_dct.
596	!!
597	!!-------------------------------------------------------------------------------------------
598	!! * Arguments
599	TYPE(SECTION),INTENT(INOUT) :: sec
600	LOGICAL ,INTENT(IN) :: ld_debug
601	INTEGER ,INTENT(IN) :: jsec ! numeric identifier of section
602
603	!! * Local variables
604	INTEGER :: jk, jseg, jclass,jl, &!loop on level/segment/classes/ice categories
605	isgnu, isgnv !
606	REAL(wp) :: zumid, zvmid, &!U/V velocity on a cell segment
607	zumid_ice, zvmid_ice, &!U/V ice velocity
608	zTnorm !transport of velocity through one cell's sides
609	REAL(wp) :: ztn, zsn, zrhoi, zrhop, zsshn, zdep !temperature/salinity/potential density/ssh/depth at u/v point
610
611	TYPE(POINT_SECTION) :: k
612	!!--------------------------------------------------------
613	!
614	IF( ld_debug )WRITE(numout,*)' Compute transport'
615
616	!---------------------------!
617	! COMPUTE TRANSPORT !
618	!---------------------------!
619	IF(sec%nb_point .NE. 0)THEN
620
621	!----------------------------------------------------------------------------------------------------
622	!Compute sign for velocities:
623	!
624	!convention:
625	! non horizontal section: direction + is toward left hand of section
626	! horizontal section: direction + is toward north of section
627	!
628	!
629	! slopeSection < 0 slopeSection > 0 slopeSection=inf slopeSection=0
630	! ---------------- ----------------- --------------- --------------
631	!
632	! isgnv=1 direction +
633	! ______ _____ ______
634	! \| //\| \| \| direction +
635	! \| isgnu=1 // \| \|isgnu=1 \|isgnu=1 /\|\
636	! \|_______ // ______\| \\ \| ---\ \|
637	! \| \| isgnv=-1 \\ \| \| ---/ direction + ____________
638	! \| \| __\\\| \|
639	! \| \| direction + \| isgnv=1
640	!
641	!----------------------------------------------------------------------------------------------------
642	isgnu = 1
643	IF( sec%slopeSection .GT. 0 ) THEN ; isgnv = -1
644	ELSE ; isgnv = 1
645	ENDIF
646	IF( sec%slopeSection .GE. 9999. ) isgnv = 1
647
648	IF( ld_debug )write(numout,*)"sec%slopeSection isgnu isgnv ",sec%slopeSection,isgnu,isgnv
649
650	!--------------------------------------!
651	! LOOP ON THE SEGMENT BETWEEN 2 NODES !
652	!--------------------------------------!
653	DO jseg=1,MAX(sec%nb_point-1,0)
654
655	!-------------------------------------------------------------------------------------------
656	! Select the appropriate coordinate for computing the velocity of the segment
657	!
658	! CASE(0) Case (2)
659	! ------- --------
660	! listPoint(jseg) listPoint(jseg+1) listPoint(jseg) F(i,j)
661	! F(i,j)----------V(i+1,j)-------F(i+1,j) \|
662	! \|
663	! \|
664	! \|
665	! Case (3) U(i,j)
666	! -------- \|
667	! \|
668	! listPoint(jseg+1) F(i,j+1) \|
669	! \| \|
670	! \| \|
671	! \| listPoint(jseg+1) F(i,j-1)
672	! \|
673	! \|
674	! U(i,j+1)
675	! \| Case(1)
676	! \| ------
677	! \|
678	! \| listPoint(jseg+1) listPoint(jseg)
679	! \| F(i-1,j)-----------V(i,j) -------f(jseg)
680	! listPoint(jseg) F(i,j)
681	!
682	!-------------------------------------------------------------------------------------------
683
684	SELECT CASE( sec%direction(jseg) )
685	CASE(0) ; k = sec%listPoint(jseg)
686	CASE(1) ; k = POINT_SECTION(sec%listPoint(jseg)%I+1,sec%listPoint(jseg)%J)
687	CASE(2) ; k = sec%listPoint(jseg)
688	CASE(3) ; k = POINT_SECTION(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J+1)
689	END SELECT
690
691	!---------------------------\|
692	! LOOP ON THE LEVEL \|
693	!---------------------------\|
694	DO jk = 1, mbkt(k%I,k%J) !Sum of the transport on the vertical
695	! ! compute temperature, salinity, insitu & potential density, ssh and depth at U/V point
696	SELECT CASE( sec%direction(jseg) )
697	CASE(0,1)
698	ztn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_tem) )
699	zsn = interp(k%I,k%J,jk,'V',tsn(:,:,:,jp_sal) )
700	zrhop = interp(k%I,k%J,jk,'V',rhop)
701	zrhoi = interp(k%I,k%J,jk,'V',rhd*rau0+rau0)
702	zsshn = 0.5( sshn(k%I,k%J) + sshn(k%I,k%J+1) ) vmask(k%I,k%J,1)
703	CASE(2,3)
704	ztn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_tem) )
705	zsn = interp(k%I,k%J,jk,'U',tsn(:,:,:,jp_sal) )
706	zrhop = interp(k%I,k%J,jk,'U',rhop)
707	zrhoi = interp(k%I,k%J,jk,'U',rhd*rau0+rau0)
708	zsshn = 0.5( sshn(k%I,k%J) + sshn(k%I+1,k%J) ) umask(k%I,k%J,1)
709	END SELECT
710	!
711	zdep= gdept_n(k%I,k%J,jk)
712
713	SELECT CASE( sec%direction(jseg) ) !compute velocity with the correct direction
714	CASE(0,1)
715	zumid=0._wp
716	zvmid=isgnvvn(k%I,k%J,jk)vmask(k%I,k%J,jk)
717	CASE(2,3)
718	zumid=isgnuun(k%I,k%J,jk)umask(k%I,k%J,jk)
719	zvmid=0._wp
720	END SELECT
721
722	!zTnorm=transport through one cell;
723	!velocity* cell's length * cell's thickness
724	zTnorm = zumide2u(k%I,k%J) e3u_n(k%I,k%J,jk) &
725	& + zvmide1v(k%I,k%J) e3v_n(k%I,k%J,jk)
726
727	!!gm THIS is WRONG no transport due to ssh in linear free surface case !!!!!
728	IF( ln_linssh ) THEN !add transport due to free surface
729	IF( jk==1 ) THEN
730	zTnorm = zTnorm + zumid* e2u(k%I,k%J) * zsshn * umask(k%I,k%J,jk) &
731	& + zvmid* e1v(k%I,k%J) * zsshn * vmask(k%I,k%J,jk)
732	ENDIF
733	ENDIF
734	!!gm end
735	!COMPUTE TRANSPORT
736
737	transports_3d(1,jsec,jseg,jk) = transports_3d(1,jsec,jseg,jk) + zTnorm
738
739	IF( sec%llstrpond ) THEN
740	transports_3d(2,jsec,jseg,jk) = transports_3d(2,jsec,jseg,jk) + zTnorm * ztn * zrhop * rcp
741	transports_3d(3,jsec,jseg,jk) = transports_3d(3,jsec,jseg,jk) + zTnorm * zsn * zrhop * 0.001
742	ENDIF
743
744	END DO !end of loop on the level
745
746	#if defined key_lim3
747
748	!ICE CASE
749	!------------
750	IF( sec%ll_ice_section )THEN
751	SELECT CASE (sec%direction(jseg))
752	CASE(0)
753	zumid_ice = 0
754	zvmid_ice = isgnv0.5(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
755	CASE(1)
756	zumid_ice = 0
757	zvmid_ice = isgnv0.5(v_ice(k%I,k%J+1)+v_ice(k%I+1,k%J+1))
758	CASE(2)
759	zvmid_ice = 0
760	zumid_ice = isgnu0.5(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
761	CASE(3)
762	zvmid_ice = 0
763	zumid_ice = isgnu0.5(u_ice(k%I+1,k%J)+u_ice(k%I+1,k%J+1))
764	END SELECT
765
766	zTnorm=zumid_icee2u(k%I,k%J)+zvmid_icee1v(k%I,k%J)
767
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	( h_i(sec%listPoint(jseg)%I,sec%listPoint(jseg)%J,jl) + &
773	h_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,mbkt(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_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), POINTER, DIMENSION(:):: zsumclasses ! 1D workspace
1011	!!-------------------------------------------------------------
1012	CALL wrk_alloc(nb_type_class , zsumclasses )
1013
1014	zsumclasses(:)=0._wp
1015	zslope = sec%slopeSection
1016
1017
1018	DO jclass=1,MAX(1,sec%nb_class-1)
1019
1020	classe = 'N '
1021	zbnd1 = 0._wp
1022	zbnd2 = 0._wp
1023	zsumclasses(1:nb_type_class)=zsumclasses(1:nb_type_class)+sec%transport(1:nb_type_class,jclass)
1024
1025
1026	!insitu density classes transports
1027	IF( ( sec%zsigi(jclass) .NE. 99._wp ) .AND. &
1028	( sec%zsigi(jclass+1) .NE. 99._wp ) )THEN
1029	classe = 'DI '
1030	zbnd1 = sec%zsigi(jclass)
1031	zbnd2 = sec%zsigi(jclass+1)
1032	ENDIF
1033	!potential density classes transports
1034	IF( ( sec%zsigp(jclass) .NE. 99._wp ) .AND. &
1035	( sec%zsigp(jclass+1) .NE. 99._wp ) )THEN
1036	classe = 'DP '
1037	zbnd1 = sec%zsigp(jclass)
1038	zbnd2 = sec%zsigp(jclass+1)
1039	ENDIF
1040	!depth classes transports
1041	IF( ( sec%zlay(jclass) .NE. 99._wp ) .AND. &
1042	( sec%zlay(jclass+1) .NE. 99._wp ) )THEN
1043	classe = 'Z '
1044	zbnd1 = sec%zlay(jclass)
1045	zbnd2 = sec%zlay(jclass+1)
1046	ENDIF
1047	!salinity classes transports
1048	IF( ( sec%zsal(jclass) .NE. 99._wp ) .AND. &
1049	( sec%zsal(jclass+1) .NE. 99._wp ) )THEN
1050	classe = 'S '
1051	zbnd1 = sec%zsal(jclass)
1052	zbnd2 = sec%zsal(jclass+1)
1053	ENDIF
1054	!temperature classes transports
1055	IF( ( sec%ztem(jclass) .NE. 99._wp ) .AND. &
1056	( sec%ztem(jclass+1) .NE. 99._wp ) ) THEN
1057	classe = 'T '
1058	zbnd1 = sec%ztem(jclass)
1059	zbnd2 = sec%ztem(jclass+1)
1060	ENDIF
1061
1062	!write volume transport per class
1063	WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, &
1064	jclass,classe,zbnd1,zbnd2,&
1065	sec%transport(1,jclass),sec%transport(2,jclass), &
1066	sec%transport(1,jclass)+sec%transport(2,jclass)
1067
1068	IF( sec%llstrpond )THEN
1069
1070	!write heat transport per class:
1071	WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope, &
1072	jclass,classe,zbnd1,zbnd2,&
1073	sec%transport(3,jclass)1.e-15,sec%transport(4,jclass)1.e-15, &
1074	( sec%transport(3,jclass)+sec%transport(4,jclass) )*1.e-15
1075	!write salt transport per class
1076	WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope, &
1077	jclass,classe,zbnd1,zbnd2,&
1078	sec%transport(5,jclass)1.e-9,sec%transport(6,jclass)1.e-9,&
1079	(sec%transport(5,jclass)+sec%transport(6,jclass))*1.e-9
1080	ENDIF
1081
1082	ENDDO
1083
1084	zbnd1 = 0._wp
1085	zbnd2 = 0._wp
1086	jclass=0
1087
1088	!write total volume transport
1089	WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope, &
1090	jclass,"total",zbnd1,zbnd2,&
1091	zsumclasses(1),zsumclasses(2),zsumclasses(1)+zsumclasses(2)
1092
1093	IF( sec%llstrpond )THEN
1094
1095	!write total heat transport
1096	WRITE(numdct_heat,119) ndastp,kt,ksec,sec%name,zslope, &
1097	jclass,"total",zbnd1,zbnd2,&
1098	zsumclasses(3)1.e-15,zsumclasses(4)1.e-15,&
1099	(zsumclasses(3)+zsumclasses(4) )*1.e-15
1100	!write total salt transport
1101	WRITE(numdct_salt,119) ndastp,kt,ksec,sec%name,zslope, &
1102	jclass,"total",zbnd1,zbnd2,&
1103	zsumclasses(5)1.e-9,zsumclasses(6)1.e-9,&
1104	(zsumclasses(5)+zsumclasses(6))*1.e-9
1105	ENDIF
1106
1107
1108	IF ( sec%ll_ice_section) THEN
1109	!write total ice volume transport
1110	WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,&
1111	jclass,"ice_vol",zbnd1,zbnd2,&
1112	sec%transport(7,1),sec%transport(8,1),&
1113	sec%transport(7,1)+sec%transport(8,1)
1114	!write total ice surface transport
1115	WRITE(numdct_vol,118) ndastp,kt,ksec,sec%name,zslope,&
1116	jclass,"ice_surf",zbnd1,zbnd2,&
1117	sec%transport(9,1),sec%transport(10,1), &
1118	sec%transport(9,1)+sec%transport(10,1)
1119	ENDIF
1120
1121	118 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3F12.4)
1122	119 FORMAT(I8,1X,I8,1X,I4,1X,A30,1X,f9.2,1X,I4,3X,A8,1X,2F12.4,5X,3E15.6)
1123
1124	CALL wrk_dealloc(nb_type_class , zsumclasses )
1125	!
1126	END SUBROUTINE dia_dct_wri
1127
1128
1129	FUNCTION interp(ki, kj, kk, cd_point, ptab)
1130	!!----------------------------------------------------------------------
1131	!!
1132	!! Purpose: compute temperature/salinity/density at U-point or V-point
1133	!! --------
1134	!!
1135	!! Method:
1136	!! ------
1137	!!
1138	!! ====> full step and partial step
1139	!!
1140	!!
1141	!! \| I \| I+1 \| Z=temperature/salinity/density at U-poinT
1142	!! \| \| \|
1143	!! ---------------------------------------- 1. Veritcal interpolation: compute zbis
1144	!! \| \| \| interpolation between ptab(I,J,K) and ptab(I,J,K+1)
1145	!! \| \| \| zbis =
1146	!! \| \| \| [ e3w(I+1,J,K)ptab(I,J,K) + ( e3w(I,J,K) - e3w(I+1,J,K) ) ptab(I,J,K-1) ]
1147	!! \| \| \| /[ e3w(I+1,J,K) + e3w(I,J,K) - e3w(I+1,J,K) ]
1148	!! \| \| \|
1149	!! \| \| \| 2. Horizontal interpolation: compute value at U/V point
1150	!!K-1 \| ptab(I,J,K-1) \| \| interpolation between zbis and ptab(I+1,J,K)
1151	!! \| . \| \|
1152	!! \| . \| \| interp = ( 0.5zet2zbis + 0.5zet1ptab(I+1,J,K) )/(0.5zet2+0.5zet1)
1153	!! \| . \| \|
1154	!! ------------------------------------------
1155	!! \| . \| \|
1156	!! \| . \| \|
1157	!! \| . \| \|
1158	!!K \| zbis.......U...ptab(I+1,J,K) \|
1159	!! \| . \| \|
1160	!! \| ptab(I,J,K) \| \|
1161	!! \| \|------------------\|
1162	!! \| \| partials \|
1163	!! \| \| steps \|
1164	!! -------------------------------------------
1165	!! <----zet1------><----zet2--------->
1166	!!
1167	!!
1168	!! ====> s-coordinate
1169	!!
1170	!! \| \| \| 1. Compute distance between T1 and U points: SQRT( zdep1^2 + (0.5 * zet1 )^2
1171	!! \| \| \| Compute distance between T2 and U points: SQRT( zdep2^2 + (0.5 * zet2 )^2
1172	!! \| \| ptab(I+1,J,K) \|
1173	!! \| \| T2 \| 2. Interpolation between T1 and T2 values at U point
1174	!! \| \| ^ \|
1175	!! \| \| \| zdep2 \|
1176	!! \| \| \| \|
1177	!! \| ^ U v \|
1178	!! \| \| \| \|
1179	!! \| \| zdep1 \| \|
1180	!! \| v \| \|
1181	!! \| T1 \| \|
1182	!! \| ptab(I,J,K) \| \|
1183	!! \| \| \|
1184	!! \| \| \|
1185	!!
1186	!! <----zet1--------><----zet2--------->
1187	!!
1188	!!----------------------------------------------------------------------
1189	!*arguments
1190	INTEGER, INTENT(IN) :: ki, kj, kk ! coordinate of point
1191	CHARACTER(len=1), INTENT(IN) :: cd_point ! type of point (U, V)
1192	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(IN) :: ptab ! variable to compute at (ki, kj, kk )
1193	REAL(wp) :: interp ! interpolated variable
1194
1195	!*local declations
1196	INTEGER :: ii1, ij1, ii2, ij2 ! local integer
1197	REAL(wp):: ze3t, ze3, zwgt1, zwgt2, zbis, zdepu ! local real
1198	REAL(wp):: zet1, zet2 ! weight for interpolation
1199	REAL(wp):: zdep1,zdep2 ! differences of depth
1200	REAL(wp):: zmsk ! mask value
1201	!!----------------------------------------------------------------------
1202
1203	IF( cd_point=='U' )THEN
1204	ii1 = ki ; ij1 = kj
1205	ii2 = ki+1 ; ij2 = kj
1206
1207	zet1=e1t(ii1,ij1)
1208	zet2=e1t(ii2,ij2)
1209	zmsk=umask(ii1,ij1,kk)
1210
1211
1212	ELSE ! cd_point=='V'
1213	ii1 = ki ; ij1 = kj
1214	ii2 = ki ; ij2 = kj+1
1215
1216	zet1=e2t(ii1,ij1)
1217	zet2=e2t(ii2,ij2)
1218	zmsk=vmask(ii1,ij1,kk)
1219
1220	ENDIF
1221
1222	IF( ln_sco )THEN ! s-coordinate case
1223
1224	zdepu = ( gdept_n(ii1,ij1,kk) + gdept_n(ii2,ij2,kk) ) * 0.5_wp
1225	zdep1 = gdept_n(ii1,ij1,kk) - zdepu
1226	zdep2 = gdept_n(ii2,ij2,kk) - zdepu
1227
1228	! weights
1229	zwgt1 = SQRT( ( 0.5 * zet1 ) * ( 0.5 * zet1 ) + ( zdep1 * zdep1 ) )
1230	zwgt2 = SQRT( ( 0.5 * zet2 ) * ( 0.5 * zet2 ) + ( zdep2 * zdep2 ) )
1231
1232	! result
1233	interp = zmsk * ( zwgt2 * ptab(ii1,ij1,kk) + zwgt1 * ptab(ii1,ij1,kk) ) / ( zwgt2 + zwgt1 )
1234
1235
1236	ELSE ! full step or partial step case
1237
1238	ze3t = e3t_n(ii2,ij2,kk) - e3t_n(ii1,ij1,kk)
1239	zwgt1 = ( e3w_n(ii2,ij2,kk) - e3w_n(ii1,ij1,kk) ) / e3w_n(ii2,ij2,kk)
1240	zwgt2 = ( e3w_n(ii1,ij1,kk) - e3w_n(ii2,ij2,kk) ) / e3w_n(ii1,ij1,kk)
1241
1242	IF(kk .NE. 1)THEN
1243
1244	IF( ze3t >= 0. )THEN
1245	! zbis
1246	zbis = ptab(ii2,ij2,kk) + zwgt1 * ( ptab(ii2,ij2,kk-1) - ptab(ii2,ij2,kk) )
1247	! result
1248	interp = zmsk * ( zet2 * ptab(ii1,ij1,kk) + zet1 * zbis )/( zet1 + zet2 )
1249	ELSE
1250	! zbis
1251	zbis = ptab(ii1,ij1,kk) + zwgt2 * ( ptab(ii1,ij1,kk-1) - ptab(ii1,ij2,kk) )
1252	! result
1253	interp = zmsk * ( zet2 * zbis + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 )
1254	ENDIF
1255
1256	ELSE
1257	interp = zmsk * ( zet2 * ptab(ii1,ij1,kk) + zet1 * ptab(ii2,ij2,kk) )/( zet1 + zet2 )
1258	ENDIF
1259
1260	ENDIF
1261	!
1262	END FUNCTION interp
1263
1264	#else
1265	!!----------------------------------------------------------------------
1266	!! Default option : Dummy module
1267	!!----------------------------------------------------------------------
1268	LOGICAL, PUBLIC, PARAMETER :: lk_diadct = .FALSE. !: diamht flag
1269	PUBLIC
1270	!! $Id$
1271	CONTAINS
1272
1273	SUBROUTINE dia_dct_init ! Dummy routine
1274	WRITE(,) 'dia_dct_init: You should not have seen this print! error?', kt
1275	END SUBROUTINE dia_dct_init
1276
1277	SUBROUTINE dia_dct( kt ) ! Dummy routine
1278	INTEGER, INTENT( in ) :: kt ! ocean time-step index
1279	WRITE(,) 'dia_dct: You should not have seen this print! error?', kt
1280	END SUBROUTINE dia_dct
1281	#endif
1282
1283	!!======================================================================
1284	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_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 8563

Download in other formats: