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

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

diadct.F90 in branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: branches/NERC/dev_r5518_NOC_MEDUSA_Stable/NEMOGCM/NEMO/OPA_SRC/DIA/diadct.F90 @ 5726

Last change on this file since 5726 was 5726, checked in by jpalmier, 9 years ago
JPALM -- 10-09-2015 -- add MEDUSA in the branch ; adapted TOP_SRC to MEDUSA ; remove some svn keywords in the branch
Property svn:executable set to ``*
File size: 59.8 KB

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

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

Download in other formats: