Context Navigation

diadct.F90 @ 9383

Last change on this file since 9383 was 9383, checked in by andmirek, 6 years ago
#2050 fixes and changes
Property svn:executable set to ``*
File size: 60.6 KB

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

Download in other formats: