Home

Context Navigation

source: palm/trunk/SOURCE/data_output_2d.f90 @ 1065

Last change on this file since 1065 was 1065, checked in by hoffmann, 11 years ago
cloud physics: rain sedimentation and turbulence effects
Property svn:keywords set to `Id`
File size: 71.6 KB

Line
1	SUBROUTINE data_output_2d( mode, av )
2
3	!--------------------------------------------------------------------------------!
4	! This file is part of PALM.
5	!
6	! PALM is free software: you can redistribute it and/or modify it under the terms
7	! of the GNU General Public License as published by the Free Software Foundation,
8	! either version 3 of the License, or (at your option) any later version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2012 Leibniz University Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	! Bugfix: Output of cross sections of ql
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: data_output_2d.f90 1065 2012-11-22 17:42:36Z hoffmann $
27	!
28	! 1053 2012-11-13 17:11:03Z hoffmann
29	! +qr, nr, qc and cross sections
30	!
31	! 1036 2012-10-22 13:43:42Z raasch
32	! code put under GPL (PALM 3.9)
33	!
34	! 1031 2012-10-19 14:35:30Z raasch
35	! netCDF4 without parallel file support implemented
36	!
37	! 1007 2012-09-19 14:30:36Z franke
38	! Bugfix: missing calculation of ql_vp added
39	!
40	! 978 2012-08-09 08:28:32Z fricke
41	! +z0h
42	!
43	! 790 2011-11-29 03:11:20Z raasch
44	! bugfix: calculation of 'pr' must depend on the particle weighting factor
45	!
46	! 771 2011-10-27 10:56:21Z heinze
47	! +lpt
48	!
49	! 759 2011-09-15 13:58:31Z raasch
50	! Splitting of parallel I/O
51	!
52	! 729 2011-05-26 10:33:34Z heinze
53	! Exchange ghost layers for p regardless of used pressure solver (except SOR).
54	!
55	! 691 2011-03-04 08:45:30Z maronga
56	! Replaced simulated_time by time_since_reference_point
57	!
58	! 673 2011-01-18 16:19:48Z suehring
59	! When using Multigrid or SOR solver an additional CALL exchange_horiz is
60	! is needed for pressure output.
61	!
62	! 667 2010-12-23 12:06:00Z suehring/gryschka
63	! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
64	! allocation of arrays local_2d and total_2d.
65	! Calls of exchange_horiz are modiefied.
66	!
67	! 622 2010-12-10 08:08:13Z raasch
68	! optional barriers included in order to speed up collective operations
69	!
70	! 493 2010-03-01 08:30:24Z raasch
71	! netCDF4 support (parallel output)
72	!
73	! 367 2009-08-25 08:35:52Z maronga
74	! simulated_time in netCDF output replaced by time_since_reference_point.
75	! Output of netCDF messages with aid of message handling routine.
76	! Bugfix: averaging along z is not allowed for 2d quantities (e.g. u* and z0)
77	! Output of messages replaced by message handling routine.
78	! Output of user defined 2D (XY) arrays at z=nzb+1 is now possible
79	! Bugfix: to_be_resorted => s_av for time-averaged scalars
80	! Calculation of shf* and qsws* added.
81	!
82	! 215 2008-11-18 09:54:31Z raasch
83	! Bugfix: no output of particle concentration and radius unless particles
84	! have been started
85	!
86	! 96 2007-06-04 08:07:41Z raasch
87	! Output of density and salinity
88	!
89	! 75 2007-03-22 09:54:05Z raasch
90	! Output of precipitation amount/rate and roughness length,
91	! 2nd+3rd argument removed from exchange horiz
92	!
93	! RCS Log replace by Id keyword, revision history cleaned up
94	!
95	! Revision 1.5 2006/08/22 13:50:29 raasch
96	! xz and yz cross sections now up to nzt+1
97	!
98	! Revision 1.2 2006/02/23 10:19:22 raasch
99	! Output of time-averaged data, output of averages along x, y, or z,
100	! output of user-defined quantities,
101	! section data are copied from local_pf to local_2d before they are output,
102	! output of particle concentration and mean radius,
103	! Former subroutine plot_2d renamed data_output_2d, pl2d.. renamed do2d..,
104	! anz renamed ngp, ebene renamed section, pl2d_.._anz renamed do2d_.._n
105	!
106	! Revision 1.1 1997/08/11 06:24:09 raasch
107	! Initial revision
108	!
109	!
110	! Description:
111	! ------------
112	! Data output of horizontal cross-sections in netCDF format or binary format
113	! compatible to old graphic software iso2d.
114	! Attention: The position of the sectional planes is still not always computed
115	! --------- correctly. (zu is used always)!
116	!------------------------------------------------------------------------------!
117
118	USE arrays_3d
119	USE averaging
120	USE cloud_parameters
121	USE control_parameters
122	USE cpulog
123	USE grid_variables
124	USE indices
125	USE interfaces
126	USE netcdf_control
127	USE particle_attributes
128	USE pegrid
129
130	IMPLICIT NONE
131
132	CHARACTER (LEN=2) :: do2d_mode, mode
133	CHARACTER (LEN=4) :: grid
134	CHARACTER (LEN=25) :: section_chr
135	CHARACTER (LEN=50) :: rtext
136	INTEGER :: av, ngp, file_id, i, if, is, iis, j, k, l, layer_xy, n, psi, &
137	s, sender, &
138	ind(4)
139	LOGICAL :: found, resorted, two_d
140	REAL :: mean_r, s_r3, s_r4
141	REAL, DIMENSION(:), ALLOCATABLE :: level_z
142	REAL, DIMENSION(:,:), ALLOCATABLE :: local_2d, local_2d_l
143	REAL, DIMENSION(:,:,:), ALLOCATABLE :: local_pf
144	#if defined( __parallel )
145	REAL, DIMENSION(:,:), ALLOCATABLE :: total_2d
146	#endif
147	REAL, DIMENSION(:,:,:), POINTER :: to_be_resorted
148
149	NAMELIST /LOCAL/ rtext
150
151	CALL cpu_log (log_point(3),'data_output_2d','start')
152
153	!
154	!-- Immediate return, if no output is requested (no respective sections
155	!-- found in parameter data_output)
156	IF ( mode == 'xy' .AND. .NOT. data_output_xy(av) ) RETURN
157	IF ( mode == 'xz' .AND. .NOT. data_output_xz(av) ) RETURN
158	IF ( mode == 'yz' .AND. .NOT. data_output_yz(av) ) RETURN
159
160	two_d = .FALSE. ! local variable to distinguish between output of pure 2D
161	! arrays and cross-sections of 3D arrays.
162
163	!
164	!-- Depending on the orientation of the cross-section, the respective output
165	!-- files have to be opened.
166	SELECT CASE ( mode )
167
168	CASE ( 'xy' )
169	s = 1
170	ALLOCATE( level_z(nzb:nzt+1), local_2d(nxlg:nxrg,nysg:nyng) )
171
172	!
173	!-- Parallel netCDF4/HDF5 output is done on all PEs, all other on PE0 only
174	IF ( netcdf_output .AND. ( myid == 0 .OR. netcdf_data_format > 4 ) ) &
175	THEN
176	CALL check_open( 101+av*10 )
177	ENDIF
178
179	IF ( data_output_2d_on_each_pe ) THEN
180	CALL check_open( 21 )
181	ELSE
182	IF ( myid == 0 ) THEN
183	IF ( iso2d_output ) CALL check_open( 21 )
184	#if defined( __parallel )
185	ALLOCATE( total_2d(-nbgp:nx+nbgp,-nbgp:ny+nbgp) )
186	#endif
187	ENDIF
188	ENDIF
189
190	CASE ( 'xz' )
191	s = 2
192	ALLOCATE( local_2d(nxlg:nxrg,nzb:nzt+1) )
193
194	!
195	!-- Parallel netCDF4/HDF5 output is done on all PEs, all other on PE0 only
196	IF ( netcdf_output .AND. ( myid == 0 .OR. netcdf_data_format > 4 ) ) &
197	THEN
198	CALL check_open( 102+av*10 )
199	ENDIF
200
201	IF ( data_output_2d_on_each_pe ) THEN
202	CALL check_open( 22 )
203	ELSE
204	IF ( myid == 0 ) THEN
205	IF ( iso2d_output ) CALL check_open( 22 )
206	#if defined( __parallel )
207	ALLOCATE( total_2d(-nbgp:nx+nbgp,nzb:nzt+1) )
208	#endif
209	ENDIF
210	ENDIF
211
212	CASE ( 'yz' )
213
214	s = 3
215	ALLOCATE( local_2d(nysg:nyng,nzb:nzt+1) )
216
217	!
218	!-- Parallel netCDF4/HDF5 output is done on all PEs, all other on PE0 only
219	IF ( netcdf_output .AND. ( myid == 0 .OR. netcdf_data_format > 4 ) ) &
220	THEN
221	CALL check_open( 103+av*10 )
222	ENDIF
223
224	IF ( data_output_2d_on_each_pe ) THEN
225	CALL check_open( 23 )
226	ELSE
227	IF ( myid == 0 ) THEN
228	IF ( iso2d_output ) CALL check_open( 23 )
229	#if defined( __parallel )
230	ALLOCATE( total_2d(-nbgp:ny+nbgp,nzb:nzt+1) )
231	#endif
232	ENDIF
233	ENDIF
234
235	CASE DEFAULT
236
237	message_string = 'unknown cross-section: ' // TRIM( mode )
238	CALL message( 'data_output_2d', 'PA0180', 1, 2, 0, 6, 0 )
239
240	END SELECT
241
242	!
243	!-- Allocate a temporary array for resorting (kji -> ijk).
244	ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb:nzt+1) )
245
246	!
247	!-- Loop of all variables to be written.
248	!-- Output dimensions chosen
249	if = 1
250	l = MAX( 2, LEN_TRIM( do2d(av,if) ) )
251	do2d_mode = do2d(av,if)(l-1:l)
252
253	DO WHILE ( do2d(av,if)(1:1) /= ' ' )
254
255	IF ( do2d_mode == mode ) THEN
256	!
257	!-- Store the array chosen on the temporary array.
258	resorted = .FALSE.
259	SELECT CASE ( TRIM( do2d(av,if) ) )
260
261	CASE ( 'e_xy', 'e_xz', 'e_yz' )
262	IF ( av == 0 ) THEN
263	to_be_resorted => e
264	ELSE
265	to_be_resorted => e_av
266	ENDIF
267	IF ( mode == 'xy' ) level_z = zu
268
269	CASE ( 'lpt_xy', 'lpt_xz', 'lpt_yz' )
270	IF ( av == 0 ) THEN
271	to_be_resorted => pt
272	ELSE
273	to_be_resorted => lpt_av
274	ENDIF
275	IF ( mode == 'xy' ) level_z = zu
276
277	CASE ( 'lwp*_xy' ) ! 2d-array
278	IF ( av == 0 ) THEN
279	DO i = nxlg, nxrg
280	DO j = nysg, nyng
281	local_pf(i,j,nzb+1) = SUM( ql(nzb:nzt,j,i) * &
282	dzw(1:nzt+1) )
283	ENDDO
284	ENDDO
285	ELSE
286	DO i = nxlg, nxrg
287	DO j = nysg, nyng
288	local_pf(i,j,nzb+1) = lwp_av(j,i)
289	ENDDO
290	ENDDO
291	ENDIF
292	resorted = .TRUE.
293	two_d = .TRUE.
294	level_z(nzb+1) = zu(nzb+1)
295
296	CASE ( 'nr_xy', 'nr_xz', 'nr_yz' )
297	IF ( av == 0 ) THEN
298	to_be_resorted => nr
299	ELSE
300	to_be_resorted => nr_av
301	ENDIF
302	IF ( mode == 'xy' ) level_z = zu
303
304	CASE ( 'p_xy', 'p_xz', 'p_yz' )
305	IF ( av == 0 ) THEN
306	IF ( psolver /= 'sor' ) CALL exchange_horiz( p, nbgp )
307	to_be_resorted => p
308	ELSE
309	IF ( psolver /= 'sor' ) CALL exchange_horiz( p_av, nbgp )
310	to_be_resorted => p_av
311	ENDIF
312	IF ( mode == 'xy' ) level_z = zu
313
314	CASE ( 'pc_xy', 'pc_xz', 'pc_yz' ) ! particle concentration
315	IF ( av == 0 ) THEN
316	IF ( simulated_time >= particle_advection_start ) THEN
317	tend = prt_count
318	CALL exchange_horiz( tend, nbgp )
319	ELSE
320	tend = 0.0
321	ENDIF
322	DO i = nxlg, nxrg
323	DO j = nysg, nyng
324	DO k = nzb, nzt+1
325	local_pf(i,j,k) = tend(k,j,i)
326	ENDDO
327	ENDDO
328	ENDDO
329	resorted = .TRUE.
330	ELSE
331	CALL exchange_horiz( pc_av, nbgp )
332	to_be_resorted => pc_av
333	ENDIF
334
335	CASE ( 'pr_xy', 'pr_xz', 'pr_yz' ) ! mean particle radius
336	IF ( av == 0 ) THEN
337	IF ( simulated_time >= particle_advection_start ) THEN
338	DO i = nxl, nxr
339	DO j = nys, nyn
340	DO k = nzb, nzt+1
341	psi = prt_start_index(k,j,i)
342	s_r3 = 0.0
343	s_r4 = 0.0
344	DO n = psi, psi+prt_count(k,j,i)-1
345	s_r3 = s_r3 + particles(n)%radius*3 &
346	particles(n)%weight_factor
347	s_r4 = s_r4 + particles(n)%radius*4 &
348	particles(n)%weight_factor
349	ENDDO
350	IF ( s_r3 /= 0.0 ) THEN
351	mean_r = s_r4 / s_r3
352	ELSE
353	mean_r = 0.0
354	ENDIF
355	tend(k,j,i) = mean_r
356	ENDDO
357	ENDDO
358	ENDDO
359	CALL exchange_horiz( tend, nbgp )
360	ELSE
361	tend = 0.0
362	END IF
363	DO i = nxlg, nxrg
364	DO j = nysg, nyng
365	DO k = nzb, nzt+1
366	local_pf(i,j,k) = tend(k,j,i)
367	ENDDO
368	ENDDO
369	ENDDO
370	resorted = .TRUE.
371	ELSE
372	CALL exchange_horiz( pr_av, nbgp )
373	to_be_resorted => pr_av
374	ENDIF
375
376	CASE ( 'pra*_xy' ) ! 2d-array / integral quantity => no av
377	CALL exchange_horiz_2d( precipitation_amount )
378	DO i = nxlg, nxrg
379	DO j = nysg, nyng
380	local_pf(i,j,nzb+1) = precipitation_amount(j,i)
381	ENDDO
382	ENDDO
383	precipitation_amount = 0.0 ! reset for next integ. interval
384	resorted = .TRUE.
385	two_d = .TRUE.
386	level_z(nzb+1) = zu(nzb+1)
387
388	CASE ( 'prr*_xy' ) ! 2d-array
389	IF ( icloud_scheme == 1 ) THEN
390	IF ( av == 0 ) THEN
391	CALL exchange_horiz_2d( precipitation_rate )
392	DO i = nxlg, nxrg
393	DO j = nysg, nyng
394	local_pf(i,j,nzb+1) = precipitation_rate(j,i)
395	ENDDO
396	ENDDO
397	ELSE
398	CALL exchange_horiz_2d( precipitation_rate_av )
399	DO i = nxlg, nxrg
400	DO j = nysg, nyng
401	local_pf(i,j,nzb+1) = precipitation_rate_av(j,i)
402	ENDDO
403	ENDDO
404	ENDIF
405	ELSE
406	IF ( av == 0 ) THEN
407	CALL exchange_horiz_2d( prr(nzb+1,:,:) )
408	DO i = nxlg, nxrg
409	DO j = nysg, nyng
410	local_pf(i,j,nzb+1) = prr(nzb+1,j,i) * hyrho(nzb+1)
411	ENDDO
412	ENDDO
413	ELSE
414	CALL exchange_horiz_2d( prr_av(nzb+1,:,:) )
415	DO i = nxlg, nxrg
416	DO j = nysg, nyng
417	local_pf(i,j,nzb+1) = prr_av(nzb+1,j,i) * hyrho(nzb+1)
418	ENDDO
419	ENDDO
420	ENDIF
421	ENDIF
422	resorted = .TRUE.
423	two_d = .TRUE.
424	level_z(nzb+1) = zu(nzb+1)
425
426	CASE ( 'prr_xy', 'prr_xz', 'prr_yz' )
427	IF ( av == 0 ) THEN
428	CALL exchange_horiz( prr, nbgp )
429	DO i = nxlg, nxrg
430	DO j = nysg, nyng
431	DO k = nzb, nzt+1
432	local_pf(i,j,k) = prr(k,j,i)
433	ENDDO
434	ENDDO
435	ENDDO
436	ELSE
437	CALL exchange_horiz( prr_av, nbgp )
438	DO i = nxlg, nxrg
439	DO j = nysg, nyng
440	DO k = nzb, nzt+1
441	local_pf(i,j,k) = prr_av(k,j,i)
442	ENDDO
443	ENDDO
444	ENDDO
445	ENDIF
446	resorted = .TRUE.
447	IF ( mode == 'xy' ) level_z = zu
448
449	CASE ( 'pt_xy', 'pt_xz', 'pt_yz' )
450	IF ( av == 0 ) THEN
451	IF ( .NOT. cloud_physics ) THEN
452	to_be_resorted => pt
453	ELSE
454	DO i = nxlg, nxrg
455	DO j = nysg, nyng
456	DO k = nzb, nzt+1
457	local_pf(i,j,k) = pt(k,j,i) + l_d_cp * &
458	pt_d_t(k) * &
459	ql(k,j,i)
460	ENDDO
461	ENDDO
462	ENDDO
463	resorted = .TRUE.
464	ENDIF
465	ELSE
466	to_be_resorted => pt_av
467	ENDIF
468	IF ( mode == 'xy' ) level_z = zu
469
470	CASE ( 'q_xy', 'q_xz', 'q_yz' )
471	IF ( av == 0 ) THEN
472	to_be_resorted => q
473	ELSE
474	to_be_resorted => q_av
475	ENDIF
476	IF ( mode == 'xy' ) level_z = zu
477
478	CASE ( 'qc_xy', 'qc_xz', 'qc_yz' )
479	IF ( av == 0 ) THEN
480	to_be_resorted => ql
481	ELSE
482	to_be_resorted => ql_av
483	ENDIF
484	IF ( mode == 'xy' ) level_z = zu
485
486	CASE ( 'ql_xy', 'ql_xz', 'ql_yz' )
487	IF ( av == 0 ) THEN
488	IF ( icloud_scheme == 0 ) THEN
489	DO i = nxlg, nxrg
490	DO j = nysg, nyng
491	DO k = nzb, nzt+1
492	local_pf(i,j,k) = ql(k,j,i) + qr(k,j,i)
493	ENDDO
494	ENDDO
495	ENDDO
496	resorted = .TRUE.
497	ELSE
498	to_be_resorted => ql
499	ENDIF
500	ELSE
501	IF ( icloud_scheme == 0 ) THEN
502	DO i = nxlg, nxrg
503	DO j = nysg, nyng
504	DO k = nzb, nzt+1
505	local_pf(i,j,k) = ql_av(k,j,i) + qr_av(k,j,i)
506	ENDDO
507	ENDDO
508	ENDDO
509	resorted = .TRUE.
510	ELSE
511	to_be_resorted => ql_av
512	ENDIF
513	ENDIF
514	IF ( mode == 'xy' ) level_z = zu
515
516	CASE ( 'ql_c_xy', 'ql_c_xz', 'ql_c_yz' )
517	IF ( av == 0 ) THEN
518	to_be_resorted => ql_c
519	ELSE
520	to_be_resorted => ql_c_av
521	ENDIF
522	IF ( mode == 'xy' ) level_z = zu
523
524	CASE ( 'ql_v_xy', 'ql_v_xz', 'ql_v_yz' )
525	IF ( av == 0 ) THEN
526	to_be_resorted => ql_v
527	ELSE
528	to_be_resorted => ql_v_av
529	ENDIF
530	IF ( mode == 'xy' ) level_z = zu
531
532	CASE ( 'ql_vp_xy', 'ql_vp_xz', 'ql_vp_yz' )
533	IF ( av == 0 ) THEN
534	IF ( simulated_time >= particle_advection_start ) THEN
535	DO i = nxl, nxr
536	DO j = nys, nyn
537	DO k = nzb, nzt+1
538	psi = prt_start_index(k,j,i)
539	DO n = psi, psi+prt_count(k,j,i)-1
540	tend(k,j,i) = tend(k,j,i) + &
541	particles(n)%weight_factor / &
542	prt_count(k,j,i)
543	ENDDO
544	ENDDO
545	ENDDO
546	ENDDO
547	CALL exchange_horiz( tend, nbgp )
548	ELSE
549	tend = 0.0
550	END IF
551	DO i = nxlg, nxrg
552	DO j = nysg, nyng
553	DO k = nzb, nzt+1
554	local_pf(i,j,k) = tend(k,j,i)
555	ENDDO
556	ENDDO
557	ENDDO
558	resorted = .TRUE.
559	ELSE
560	CALL exchange_horiz( ql_vp_av, nbgp )
561	to_be_resorted => ql_vp
562	ENDIF
563	IF ( mode == 'xy' ) level_z = zu
564
565	CASE ( 'qr_xy', 'qr_xz', 'qr_yz' )
566	IF ( av == 0 ) THEN
567	to_be_resorted => qr
568	ELSE
569	to_be_resorted => qr_av
570	ENDIF
571	IF ( mode == 'xy' ) level_z = zu
572
573	CASE ( 'qsws*_xy' ) ! 2d-array
574	IF ( av == 0 ) THEN
575	DO i = nxlg, nxrg
576	DO j = nysg, nyng
577	local_pf(i,j,nzb+1) = qsws(j,i)
578	ENDDO
579	ENDDO
580	ELSE
581	DO i = nxlg, nxrg
582	DO j = nysg, nyng
583	local_pf(i,j,nzb+1) = qsws_av(j,i)
584	ENDDO
585	ENDDO
586	ENDIF
587	resorted = .TRUE.
588	two_d = .TRUE.
589	level_z(nzb+1) = zu(nzb+1)
590
591	CASE ( 'qv_xy', 'qv_xz', 'qv_yz' )
592	IF ( av == 0 ) THEN
593	DO i = nxlg, nxrg
594	DO j = nysg, nyng
595	DO k = nzb, nzt+1
596	local_pf(i,j,k) = q(k,j,i) - ql(k,j,i)
597	ENDDO
598	ENDDO
599	ENDDO
600	resorted = .TRUE.
601	ELSE
602	to_be_resorted => qv_av
603	ENDIF
604	IF ( mode == 'xy' ) level_z = zu
605
606	CASE ( 'rho_xy', 'rho_xz', 'rho_yz' )
607	IF ( av == 0 ) THEN
608	to_be_resorted => rho
609	ELSE
610	to_be_resorted => rho_av
611	ENDIF
612
613	CASE ( 's_xy', 's_xz', 's_yz' )
614	IF ( av == 0 ) THEN
615	to_be_resorted => q
616	ELSE
617	to_be_resorted => s_av
618	ENDIF
619
620	CASE ( 'sa_xy', 'sa_xz', 'sa_yz' )
621	IF ( av == 0 ) THEN
622	to_be_resorted => sa
623	ELSE
624	to_be_resorted => sa_av
625	ENDIF
626
627	CASE ( 'shf*_xy' ) ! 2d-array
628	IF ( av == 0 ) THEN
629	DO i = nxlg, nxrg
630	DO j = nysg, nyng
631	local_pf(i,j,nzb+1) = shf(j,i)
632	ENDDO
633	ENDDO
634	ELSE
635	DO i = nxlg, nxrg
636	DO j = nysg, nyng
637	local_pf(i,j,nzb+1) = shf_av(j,i)
638	ENDDO
639	ENDDO
640	ENDIF
641	resorted = .TRUE.
642	two_d = .TRUE.
643	level_z(nzb+1) = zu(nzb+1)
644
645	CASE ( 't*_xy' ) ! 2d-array
646	IF ( av == 0 ) THEN
647	DO i = nxlg, nxrg
648	DO j = nysg, nyng
649	local_pf(i,j,nzb+1) = ts(j,i)
650	ENDDO
651	ENDDO
652	ELSE
653	DO i = nxlg, nxrg
654	DO j = nysg, nyng
655	local_pf(i,j,nzb+1) = ts_av(j,i)
656	ENDDO
657	ENDDO
658	ENDIF
659	resorted = .TRUE.
660	two_d = .TRUE.
661	level_z(nzb+1) = zu(nzb+1)
662
663	CASE ( 'u_xy', 'u_xz', 'u_yz' )
664	IF ( av == 0 ) THEN
665	to_be_resorted => u
666	ELSE
667	to_be_resorted => u_av
668	ENDIF
669	IF ( mode == 'xy' ) level_z = zu
670	!
671	!-- Substitute the values generated by "mirror" boundary condition
672	!-- at the bottom boundary by the real surface values.
673	IF ( do2d(av,if) == 'u_xz' .OR. do2d(av,if) == 'u_yz' ) THEN
674	IF ( ibc_uv_b == 0 ) local_pf(:,:,nzb) = 0.0
675	ENDIF
676
677	CASE ( 'u*_xy' ) ! 2d-array
678	IF ( av == 0 ) THEN
679	DO i = nxlg, nxrg
680	DO j = nysg, nyng
681	local_pf(i,j,nzb+1) = us(j,i)
682	ENDDO
683	ENDDO
684	ELSE
685	DO i = nxlg, nxrg
686	DO j = nysg, nyng
687	local_pf(i,j,nzb+1) = us_av(j,i)
688	ENDDO
689	ENDDO
690	ENDIF
691	resorted = .TRUE.
692	two_d = .TRUE.
693	level_z(nzb+1) = zu(nzb+1)
694
695	CASE ( 'v_xy', 'v_xz', 'v_yz' )
696	IF ( av == 0 ) THEN
697	to_be_resorted => v
698	ELSE
699	to_be_resorted => v_av
700	ENDIF
701	IF ( mode == 'xy' ) level_z = zu
702	!
703	!-- Substitute the values generated by "mirror" boundary condition
704	!-- at the bottom boundary by the real surface values.
705	IF ( do2d(av,if) == 'v_xz' .OR. do2d(av,if) == 'v_yz' ) THEN
706	IF ( ibc_uv_b == 0 ) local_pf(:,:,nzb) = 0.0
707	ENDIF
708
709	CASE ( 'vpt_xy', 'vpt_xz', 'vpt_yz' )
710	IF ( av == 0 ) THEN
711	to_be_resorted => vpt
712	ELSE
713	to_be_resorted => vpt_av
714	ENDIF
715	IF ( mode == 'xy' ) level_z = zu
716
717	CASE ( 'w_xy', 'w_xz', 'w_yz' )
718	IF ( av == 0 ) THEN
719	to_be_resorted => w
720	ELSE
721	to_be_resorted => w_av
722	ENDIF
723	IF ( mode == 'xy' ) level_z = zw
724
725	CASE ( 'z0*_xy' ) ! 2d-array
726	IF ( av == 0 ) THEN
727	DO i = nxlg, nxrg
728	DO j = nysg, nyng
729	local_pf(i,j,nzb+1) = z0(j,i)
730	ENDDO
731	ENDDO
732	ELSE
733	DO i = nxlg, nxrg
734	DO j = nysg, nyng
735	local_pf(i,j,nzb+1) = z0_av(j,i)
736	ENDDO
737	ENDDO
738	ENDIF
739	resorted = .TRUE.
740	two_d = .TRUE.
741	level_z(nzb+1) = zu(nzb+1)
742
743	CASE ( 'z0h*_xy' ) ! 2d-array
744	IF ( av == 0 ) THEN
745	DO i = nxlg, nxrg
746	DO j = nysg, nyng
747	local_pf(i,j,nzb+1) = z0h(j,i)
748	ENDDO
749	ENDDO
750	ELSE
751	DO i = nxlg, nxrg
752	DO j = nysg, nyng
753	local_pf(i,j,nzb+1) = z0h_av(j,i)
754	ENDDO
755	ENDDO
756	ENDIF
757	resorted = .TRUE.
758	two_d = .TRUE.
759	level_z(nzb+1) = zu(nzb+1)
760
761	CASE DEFAULT
762	!
763	!-- User defined quantity
764	CALL user_data_output_2d( av, do2d(av,if), found, grid, &
765	local_pf, two_d )
766	resorted = .TRUE.
767
768	IF ( grid == 'zu' ) THEN
769	IF ( mode == 'xy' ) level_z = zu
770	ELSEIF ( grid == 'zw' ) THEN
771	IF ( mode == 'xy' ) level_z = zw
772	ELSEIF ( grid == 'zu1' ) THEN
773	IF ( mode == 'xy' ) level_z(nzb+1) = zu(nzb+1)
774	ENDIF
775
776	IF ( .NOT. found ) THEN
777	message_string = 'no output provided for: ' // &
778	TRIM( do2d(av,if) )
779	CALL message( 'data_output_2d', 'PA0181', 0, 0, 0, 6, 0 )
780	ENDIF
781
782	END SELECT
783
784	!
785	!-- Resort the array to be output, if not done above
786	IF ( .NOT. resorted ) THEN
787	DO i = nxlg, nxrg
788	DO j = nysg, nyng
789	DO k = nzb, nzt+1
790	local_pf(i,j,k) = to_be_resorted(k,j,i)
791	ENDDO
792	ENDDO
793	ENDDO
794	ENDIF
795
796	!
797	!-- Output of the individual cross-sections, depending on the cross-
798	!-- section mode chosen.
799	is = 1
800	loop1: DO WHILE ( section(is,s) /= -9999 .OR. two_d )
801
802	SELECT CASE ( mode )
803
804	CASE ( 'xy' )
805	!
806	!-- Determine the cross section index
807	IF ( two_d ) THEN
808	layer_xy = nzb+1
809	ELSE
810	layer_xy = section(is,s)
811	ENDIF
812
813	!
814	!-- Update the netCDF xy cross section time axis
815	IF ( myid == 0 .OR. netcdf_data_format > 4 ) THEN
816	IF ( simulated_time /= do2d_xy_last_time(av) ) THEN
817	do2d_xy_time_count(av) = do2d_xy_time_count(av) + 1
818	do2d_xy_last_time(av) = simulated_time
819	IF ( ( .NOT. data_output_2d_on_each_pe .AND. &
820	netcdf_output ) .OR. netcdf_data_format > 4 ) &
821	THEN
822	#if defined( __netcdf )
823	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
824	id_var_time_xy(av), &
825	(/ time_since_reference_point /), &
826	start = (/ do2d_xy_time_count(av) /), &
827	count = (/ 1 /) )
828	CALL handle_netcdf_error( 'data_output_2d', 53 )
829	#endif
830	ENDIF
831	ENDIF
832	ENDIF
833	!
834	!-- If required, carry out averaging along z
835	IF ( section(is,s) == -1 .AND. .NOT. two_d ) THEN
836
837	local_2d = 0.0
838	!
839	!-- Carry out the averaging (all data are on the PE)
840	DO k = nzb, nzt+1
841	DO j = nysg, nyng
842	DO i = nxlg, nxrg
843	local_2d(i,j) = local_2d(i,j) + local_pf(i,j,k)
844	ENDDO
845	ENDDO
846	ENDDO
847
848	local_2d = local_2d / ( nzt -nzb + 2.0)
849
850	ELSE
851	!
852	!-- Just store the respective section on the local array
853	local_2d = local_pf(:,:,layer_xy)
854
855	ENDIF
856
857	#if defined( __parallel )
858	IF ( netcdf_output .AND. netcdf_data_format > 4 ) THEN
859	!
860	!-- Parallel output in netCDF4/HDF5 format.
861	!-- Do not output redundant ghost point data except for the
862	!-- boundaries of the total domain.
863	IF ( two_d ) THEN
864	iis = 1
865	ELSE
866	iis = is
867	ENDIF
868
869	#if defined( __netcdf )
870	IF ( nxr == nx .AND. nyn /= ny ) THEN
871	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
872	id_var_do2d(av,if), &
873	local_2d(nxl:nxr+1,nys:nyn), &
874	start = (/ nxl+1, nys+1, iis, &
875	do2d_xy_time_count(av) /), &
876	count = (/ nxr-nxl+2, &
877	nyn-nys+1, 1, 1 /) )
878	ELSEIF ( nxr /= nx .AND. nyn == ny ) THEN
879	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
880	id_var_do2d(av,if), &
881	local_2d(nxl:nxr,nys:nyn+1), &
882	start = (/ nxl+1, nys+1, iis, &
883	do2d_xy_time_count(av) /), &
884	count = (/ nxr-nxl+1, &
885	nyn-nys+2, 1, 1 /) )
886	ELSEIF ( nxr == nx .AND. nyn == ny ) THEN
887	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
888	id_var_do2d(av,if), &
889	local_2d(nxl:nxr+1,nys:nyn+1),&
890	start = (/ nxl+1, nys+1, iis, &
891	do2d_xy_time_count(av) /), &
892	count = (/ nxr-nxl+2, &
893	nyn-nys+2, 1, 1 /) )
894	ELSE
895	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
896	id_var_do2d(av,if), &
897	local_2d(nxl:nxr,nys:nyn), &
898	start = (/ nxl+1, nys+1, iis, &
899	do2d_xy_time_count(av) /), &
900	count = (/ nxr-nxl+1, &
901	nyn-nys+1, 1, 1 /) )
902	ENDIF
903
904	CALL handle_netcdf_error( 'data_output_2d', 55 )
905	#endif
906	ELSE
907
908	IF ( data_output_2d_on_each_pe ) THEN
909	!
910	!-- Output of partial arrays on each PE
911	#if defined( __netcdf )
912	IF ( netcdf_output .AND. myid == 0 ) THEN
913	WRITE ( 21 ) time_since_reference_point, &
914	do2d_xy_time_count(av), av
915	ENDIF
916	#endif
917	DO i = 0, io_blocks-1
918	IF ( i == io_group ) THEN
919	WRITE ( 21 ) nxlg, nxrg, nysg, nyng
920	WRITE ( 21 ) local_2d
921	ENDIF
922	#if defined( __parallel )
923	CALL MPI_BARRIER( comm2d, ierr )
924	#endif
925	ENDDO
926
927	ELSE
928	!
929	!-- PE0 receives partial arrays from all processors and
930	!-- then outputs them. Here a barrier has to be set,
931	!-- because otherwise "-MPI- FATAL: Remote protocol queue
932	!-- full" may occur.
933	CALL MPI_BARRIER( comm2d, ierr )
934
935	ngp = ( nxrg-nxlg+1 ) * ( nyng-nysg+1 )
936	IF ( myid == 0 ) THEN
937	!
938	!-- Local array can be relocated directly.
939	total_2d(nxlg:nxrg,nysg:nyng) = local_2d
940	!
941	!-- Receive data from all other PEs.
942	DO n = 1, numprocs-1
943	!
944	!-- Receive index limits first, then array.
945	!-- Index limits are received in arbitrary order from
946	!-- the PEs.
947	CALL MPI_RECV( ind(1), 4, MPI_INTEGER, &
948	MPI_ANY_SOURCE, 0, comm2d, &
949	status, ierr )
950	sender = status(MPI_SOURCE)
951	DEALLOCATE( local_2d )
952	ALLOCATE( local_2d(ind(1):ind(2),ind(3):ind(4)) )
953	CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, &
954	MPI_REAL, sender, 1, comm2d, &
955	status, ierr )
956	total_2d(ind(1):ind(2),ind(3):ind(4)) = local_2d
957	ENDDO
958	!
959	!-- Output of the total cross-section.
960	IF ( iso2d_output ) THEN
961	WRITE (21) total_2d(-nbgp:nx+nbgp,-nbgp:ny+nbgp)
962	ENDIF
963	!
964	!-- Relocate the local array for the next loop increment
965	DEALLOCATE( local_2d )
966	ALLOCATE( local_2d(nxlg:nxrg,nysg:nyng) )
967
968	#if defined( __netcdf )
969	IF ( netcdf_output ) THEN
970	IF ( two_d ) THEN
971	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
972	id_var_do2d(av,if), &
973	total_2d(0:nx+1,0:ny+1), &
974	start = (/ 1, 1, 1, do2d_xy_time_count(av) /), &
975	count = (/ nx+2, ny+2, 1, 1 /) )
976	ELSE
977	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
978	id_var_do2d(av,if), &
979	total_2d(0:nx+1,0:ny+1), &
980	start = (/ 1, 1, is, do2d_xy_time_count(av) /), &
981	count = (/ nx+2, ny+2, 1, 1 /) )
982	ENDIF
983	CALL handle_netcdf_error( 'data_output_2d', 54 )
984	ENDIF
985	#endif
986
987	ELSE
988	!
989	!-- First send the local index limits to PE0
990	ind(1) = nxlg; ind(2) = nxrg
991	ind(3) = nysg; ind(4) = nyng
992	CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, &
993	comm2d, ierr )
994	!
995	!-- Send data to PE0
996	CALL MPI_SEND( local_2d(nxlg,nysg), ngp, &
997	MPI_REAL, 0, 1, comm2d, ierr )
998	ENDIF
999	!
1000	!-- A barrier has to be set, because otherwise some PEs may
1001	!-- proceed too fast so that PE0 may receive wrong data on
1002	!-- tag 0
1003	CALL MPI_BARRIER( comm2d, ierr )
1004	ENDIF
1005
1006	ENDIF
1007	#else
1008	IF ( iso2d_output ) THEN
1009	WRITE (21) local_2d(nxl:nxr+1,nys:nyn+1)
1010	ENDIF
1011	#if defined( __netcdf )
1012	IF ( netcdf_output ) THEN
1013	IF ( two_d ) THEN
1014	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1015	id_var_do2d(av,if), &
1016	local_2d(nxl:nxr+1,nys:nyn+1), &
1017	start = (/ 1, 1, 1, do2d_xy_time_count(av) /), &
1018	count = (/ nx+2, ny+2, 1, 1 /) )
1019	ELSE
1020	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
1021	id_var_do2d(av,if), &
1022	local_2d(nxl:nxr+1,nys:nyn+1), &
1023	start = (/ 1, 1, is, do2d_xy_time_count(av) /), &
1024	count = (/ nx+2, ny+2, 1, 1 /) )
1025	ENDIF
1026	CALL handle_netcdf_error( 'data_output_2d', 447 )
1027	ENDIF
1028	#endif
1029	#endif
1030	do2d_xy_n = do2d_xy_n + 1
1031	!
1032	!-- Write LOCAL parameter set for ISO2D.
1033	IF ( myid == 0 .AND. iso2d_output ) THEN
1034	IF ( section(is,s) /= -1 ) THEN
1035	WRITE ( section_chr, '(''z = '',F7.2,'' m (GP '',I3, &
1036	&'')'')' &
1037	) level_z(layer_xy), layer_xy
1038	ELSE
1039	section_chr = 'averaged along z'
1040	ENDIF
1041	IF ( av == 0 ) THEN
1042	rtext = TRIM( do2d(av,if) ) // ' t = ' // &
1043	TRIM( simulated_time_chr ) // ' ' // &
1044	TRIM( section_chr )
1045	ELSE
1046	rtext = TRIM( do2d(av,if) ) // ' averaged t = ' // &
1047	TRIM( simulated_time_chr ) // ' ' // &
1048	TRIM( section_chr )
1049	ENDIF
1050	WRITE (27,LOCAL)
1051	ENDIF
1052	!
1053	!-- For 2D-arrays (e.g. u*) only one cross-section is available.
1054	!-- Hence exit loop of output levels.
1055	IF ( two_d ) THEN
1056	two_d = .FALSE.
1057	EXIT loop1
1058	ENDIF
1059
1060	CASE ( 'xz' )
1061	!
1062	!-- Update the netCDF xz cross section time axis
1063	IF ( myid == 0 .OR. netcdf_data_format > 4 ) THEN
1064
1065	IF ( simulated_time /= do2d_xz_last_time(av) ) THEN
1066	do2d_xz_time_count(av) = do2d_xz_time_count(av) + 1
1067	do2d_xz_last_time(av) = simulated_time
1068	IF ( ( .NOT. data_output_2d_on_each_pe .AND. &
1069	netcdf_output ) .OR. netcdf_data_format > 4 ) &
1070	THEN
1071	#if defined( __netcdf )
1072	nc_stat = NF90_PUT_VAR( id_set_xz(av), &
1073	id_var_time_xz(av), &
1074	(/ time_since_reference_point /), &
1075	start = (/ do2d_xz_time_count(av) /), &
1076	count = (/ 1 /) )
1077	CALL handle_netcdf_error( 'data_output_2d', 56 )
1078	#endif
1079	ENDIF
1080	ENDIF
1081
1082	ENDIF
1083
1084	!
1085	!-- If required, carry out averaging along y
1086	IF ( section(is,s) == -1 ) THEN
1087
1088	ALLOCATE( local_2d_l(nxlg:nxrg,nzb:nzt+1) )
1089	local_2d_l = 0.0
1090	ngp = ( nxrg-nxlg+1 ) * ( nzt-nzb+2 )
1091	!
1092	!-- First local averaging on the PE
1093	DO k = nzb, nzt+1
1094	DO j = nys, nyn
1095	DO i = nxlg, nxrg
1096	local_2d_l(i,k) = local_2d_l(i,k) + &
1097	local_pf(i,j,k)
1098	ENDDO
1099	ENDDO
1100	ENDDO
1101	#if defined( __parallel )
1102	!
1103	!-- Now do the averaging over all PEs along y
1104	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1105	CALL MPI_ALLREDUCE( local_2d_l(nxlg,nzb), &
1106	local_2d(nxlg,nzb), ngp, MPI_REAL, &
1107	MPI_SUM, comm1dy, ierr )
1108	#else
1109	local_2d = local_2d_l
1110	#endif
1111	local_2d = local_2d / ( ny + 1.0 )
1112
1113	DEALLOCATE( local_2d_l )
1114
1115	ELSE
1116	!
1117	!-- Just store the respective section on the local array
1118	!-- (but only if it is available on this PE!)
1119	IF ( section(is,s) >= nys .AND. section(is,s) <= nyn ) &
1120	THEN
1121	local_2d = local_pf(:,section(is,s),nzb:nzt+1)
1122	ENDIF
1123
1124	ENDIF
1125
1126	#if defined( __parallel )
1127	IF ( netcdf_output .AND. netcdf_data_format > 4 ) THEN
1128	!
1129	!-- ATTENTION: The following lines are a workaround, because
1130	!-- independet output does not work with the
1131	!-- current netCDF4 installation. Therefore, data
1132	!-- are transferred from PEs having the cross
1133	!-- sections to other PEs along y having no cross
1134	!-- section data. Some of these data are the
1135	!-- output.
1136	!-- BEGIN WORKAROUND---------------------------------------
1137	IF ( npey /= 1 .AND. section(is,s) /= -1) THEN
1138	ALLOCATE( local_2d_l(nxlg:nxrg,nzb:nzt+1) )
1139	local_2d_l = 0.0
1140	IF ( section(is,s) >= nys .AND. section(is,s) <= nyn )&
1141	THEN
1142	local_2d_l = local_2d
1143	ENDIF
1144	#if defined( __parallel )
1145	!
1146	!-- Distribute data over all PEs along y
1147	ngp = ( nxrg-nxlg+1 ) * ( nzt-nzb+2 )
1148	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1149	CALL MPI_ALLREDUCE( local_2d_l(nxlg,nzb), &
1150	local_2d(nxlg,nzb), ngp, &
1151	MPI_REAL, MPI_SUM, comm1dy, ierr )
1152	#else
1153	local_2d = local_2d_l
1154	#endif
1155	DEALLOCATE( local_2d_l )
1156	ENDIF
1157	!-- END WORKAROUND-----------------------------------------
1158
1159	!
1160	!-- Output in netCDF4/HDF5 format.
1161	!-- Output only on those PEs where the respective cross
1162	!-- sections reside. Cross sections averaged along y are
1163	!-- output on the respective first PE along y (myidy=0).
1164	IF ( ( section(is,s) >= nys .AND. &
1165	section(is,s) <= nyn ) .OR. &
1166	( section(is,s) == -1 .AND. myidy == 0 ) ) THEN
1167	!
1168	!-- Do not output redundant ghost point data except for the
1169	!-- boundaries of the total domain.
1170	#if defined( __netcdf )
1171	IF ( nxr == nx ) THEN
1172	nc_stat = NF90_PUT_VAR( id_set_xz(av), &
1173	id_var_do2d(av,if), &
1174	local_2d(nxl:nxr+1,nzb:nzt+1), &
1175	start = (/ nxl+1, is, 1, &
1176	do2d_xz_time_count(av) /), &
1177	count = (/ nxr-nxl+2, 1, &
1178	nzt+2, 1 /) )
1179	ELSE
1180	nc_stat = NF90_PUT_VAR( id_set_xz(av), &
1181	id_var_do2d(av,if), &
1182	local_2d(nxl:nxr,nzb:nzt+1), &
1183	start = (/ nxl+1, is, 1, &
1184	do2d_xz_time_count(av) /), &
1185	count = (/ nxr-nxl+1, 1, &
1186	nzt+2, 1 /) )
1187	ENDIF
1188
1189	CALL handle_netcdf_error( 'data_output_2d', 57 )
1190
1191	ELSE
1192	!
1193	!-- Output on other PEs. Only one point is output!!
1194	!-- ATTENTION: This is a workaround (see above)!!
1195	IF ( npey /= 1 ) THEN
1196	nc_stat = NF90_PUT_VAR( id_set_xz(av), &
1197	id_var_do2d(av,if), &
1198	local_2d(nxl:nxl,nzb:nzb), &
1199	start = (/ nxl+1, is, 1, &
1200	do2d_xz_time_count(av) /), &
1201	count = (/ 1, 1, 1, 1 /) )
1202	CALL handle_netcdf_error( 'data_output_2d', 451 )
1203	ENDIF
1204	#endif
1205	ENDIF
1206
1207	ELSE
1208
1209	IF ( data_output_2d_on_each_pe ) THEN
1210	!
1211	!-- Output of partial arrays on each PE. If the cross
1212	!-- section does not reside on the PE, output special
1213	!-- index values.
1214	#if defined( __netcdf )
1215	IF ( netcdf_output .AND. myid == 0 ) THEN
1216	WRITE ( 22 ) time_since_reference_point, &
1217	do2d_xz_time_count(av), av
1218	ENDIF
1219	#endif
1220	DO i = 0, io_blocks-1
1221	IF ( i == io_group ) THEN
1222	IF ( ( section(is,s) >= nys .AND. &
1223	section(is,s) <= nyn ) .OR. &
1224	( section(is,s) == -1 .AND. &
1225	nys-1 == -1 ) ) &
1226	THEN
1227	WRITE (22) nxlg, nxrg, nzb, nzt+1
1228	WRITE (22) local_2d
1229	ELSE
1230	WRITE (22) -1, -1, -1, -1
1231	ENDIF
1232	ENDIF
1233	#if defined( __parallel )
1234	CALL MPI_BARRIER( comm2d, ierr )
1235	#endif
1236	ENDDO
1237
1238	ELSE
1239	!
1240	!-- PE0 receives partial arrays from all processors of the
1241	!-- respective cross section and outputs them. Here a
1242	!-- barrier has to be set, because otherwise
1243	!-- "-MPI- FATAL: Remote protocol queue full" may occur.
1244	CALL MPI_BARRIER( comm2d, ierr )
1245
1246	ngp = ( nxrg-nxlg+1 ) * ( nzt-nzb+2 )
1247	IF ( myid == 0 ) THEN
1248	!
1249	!-- Local array can be relocated directly.
1250	IF ( ( section(is,s) >= nys .AND. &
1251	section(is,s) <= nyn ) .OR. &
1252	( section(is,s) == -1 .AND. nys-1 == -1 ) ) &
1253	THEN
1254	total_2d(nxlg:nxrg,nzb:nzt+1) = local_2d
1255	ENDIF
1256	!
1257	!-- Receive data from all other PEs.
1258	DO n = 1, numprocs-1
1259	!
1260	!-- Receive index limits first, then array.
1261	!-- Index limits are received in arbitrary order from
1262	!-- the PEs.
1263	CALL MPI_RECV( ind(1), 4, MPI_INTEGER, &
1264	MPI_ANY_SOURCE, 0, comm2d, &
1265	status, ierr )
1266	!
1267	!-- Not all PEs have data for XZ-cross-section.
1268	IF ( ind(1) /= -9999 ) THEN
1269	sender = status(MPI_SOURCE)
1270	DEALLOCATE( local_2d )
1271	ALLOCATE( local_2d(ind(1):ind(2), &
1272	ind(3):ind(4)) )
1273	CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, &
1274	MPI_REAL, sender, 1, comm2d, &
1275	status, ierr )
1276	total_2d(ind(1):ind(2),ind(3):ind(4)) = &
1277	local_2d
1278	ENDIF
1279	ENDDO
1280	!
1281	!-- Output of the total cross-section.
1282	IF ( iso2d_output ) THEN
1283	WRITE (22) total_2d(-nbgp:nx+nbgp,nzb:nzt+1)
1284	ENDIF
1285	!
1286	!-- Relocate the local array for the next loop increment
1287	DEALLOCATE( local_2d )
1288	ALLOCATE( local_2d(nxlg:nxrg,nzb:nzt+1) )
1289
1290	#if defined( __netcdf )
1291	IF ( netcdf_output ) THEN
1292	nc_stat = NF90_PUT_VAR( id_set_xz(av), &
1293	id_var_do2d(av,if), &
1294	total_2d(0:nx+1,nzb:nzt+1),&
1295	start = (/ 1, is, 1, do2d_xz_time_count(av) /), &
1296	count = (/ nx+2, 1, nz+2, 1 /) )
1297	CALL handle_netcdf_error( 'data_output_2d', 58 )
1298	ENDIF
1299	#endif
1300
1301	ELSE
1302	!
1303	!-- If the cross section resides on the PE, send the
1304	!-- local index limits, otherwise send -9999 to PE0.
1305	IF ( ( section(is,s) >= nys .AND. &
1306	section(is,s) <= nyn ) .OR. &
1307	( section(is,s) == -1 .AND. nys-1 == -1 ) ) &
1308	THEN
1309	ind(1) = nxlg; ind(2) = nxrg
1310	ind(3) = nzb; ind(4) = nzt+1
1311	ELSE
1312	ind(1) = -9999; ind(2) = -9999
1313	ind(3) = -9999; ind(4) = -9999
1314	ENDIF
1315	CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, &
1316	comm2d, ierr )
1317	!
1318	!-- If applicable, send data to PE0.
1319	IF ( ind(1) /= -9999 ) THEN
1320	CALL MPI_SEND( local_2d(nxlg,nzb), ngp, &
1321	MPI_REAL, 0, 1, comm2d, ierr )
1322	ENDIF
1323	ENDIF
1324	!
1325	!-- A barrier has to be set, because otherwise some PEs may
1326	!-- proceed too fast so that PE0 may receive wrong data on
1327	!-- tag 0
1328	CALL MPI_BARRIER( comm2d, ierr )
1329	ENDIF
1330
1331	ENDIF
1332	#else
1333	IF ( iso2d_output ) THEN
1334	WRITE (22) local_2d(nxl:nxr+1,nzb:nzt+1)
1335	ENDIF
1336	#if defined( __netcdf )
1337	IF ( netcdf_output ) THEN
1338	nc_stat = NF90_PUT_VAR( id_set_xz(av), &
1339	id_var_do2d(av,if), &
1340	local_2d(nxl:nxr+1,nzb:nzt+1), &
1341	start = (/ 1, is, 1, do2d_xz_time_count(av) /), &
1342	count = (/ nx+2, 1, nz+2, 1 /) )
1343	CALL handle_netcdf_error( 'data_output_2d', 451 )
1344	ENDIF
1345	#endif
1346	#endif
1347	do2d_xz_n = do2d_xz_n + 1
1348	!
1349	!-- Write LOCAL-parameter set for ISO2D.
1350	IF ( myid == 0 .AND. iso2d_output ) THEN
1351	IF ( section(is,s) /= -1 ) THEN
1352	WRITE ( section_chr, '(''y = '',F8.2,'' m (GP '',I3, &
1353	&'')'')' &
1354	) section(is,s)*dy, section(is,s)
1355	ELSE
1356	section_chr = 'averaged along y'
1357	ENDIF
1358	IF ( av == 0 ) THEN
1359	rtext = TRIM( do2d(av,if) ) // ' t = ' // &
1360	TRIM( simulated_time_chr ) // ' ' // &
1361	TRIM( section_chr )
1362	ELSE
1363	rtext = TRIM( do2d(av,if) ) // ' averaged t = ' // &
1364	TRIM( simulated_time_chr ) // ' ' // &
1365	TRIM( section_chr )
1366	ENDIF
1367	WRITE (28,LOCAL)
1368	ENDIF
1369
1370	CASE ( 'yz' )
1371	!
1372	!-- Update the netCDF yz cross section time axis
1373	IF ( myid == 0 .OR. netcdf_data_format > 4 ) THEN
1374
1375	IF ( simulated_time /= do2d_yz_last_time(av) ) THEN
1376	do2d_yz_time_count(av) = do2d_yz_time_count(av) + 1
1377	do2d_yz_last_time(av) = simulated_time
1378	IF ( ( .NOT. data_output_2d_on_each_pe .AND. &
1379	netcdf_output ) .OR. netcdf_data_format > 4 ) &
1380	THEN
1381	#if defined( __netcdf )
1382	nc_stat = NF90_PUT_VAR( id_set_yz(av), &
1383	id_var_time_yz(av), &
1384	(/ time_since_reference_point /), &
1385	start = (/ do2d_yz_time_count(av) /), &
1386	count = (/ 1 /) )
1387	CALL handle_netcdf_error( 'data_output_2d', 59 )
1388	#endif
1389	ENDIF
1390	ENDIF
1391
1392	ENDIF
1393	!
1394	!-- If required, carry out averaging along x
1395	IF ( section(is,s) == -1 ) THEN
1396
1397	ALLOCATE( local_2d_l(nysg:nyng,nzb:nzt+1) )
1398	local_2d_l = 0.0
1399	ngp = ( nyng-nysg+1 ) * ( nzt-nzb+2 )
1400	!
1401	!-- First local averaging on the PE
1402	DO k = nzb, nzt+1
1403	DO j = nysg, nyng
1404	DO i = nxl, nxr
1405	local_2d_l(j,k) = local_2d_l(j,k) + &
1406	local_pf(i,j,k)
1407	ENDDO
1408	ENDDO
1409	ENDDO
1410	#if defined( __parallel )
1411	!
1412	!-- Now do the averaging over all PEs along x
1413	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1414	CALL MPI_ALLREDUCE( local_2d_l(nysg,nzb), &
1415	local_2d(nysg,nzb), ngp, MPI_REAL, &
1416	MPI_SUM, comm1dx, ierr )
1417	#else
1418	local_2d = local_2d_l
1419	#endif
1420	local_2d = local_2d / ( nx + 1.0 )
1421
1422	DEALLOCATE( local_2d_l )
1423
1424	ELSE
1425	!
1426	!-- Just store the respective section on the local array
1427	!-- (but only if it is available on this PE!)
1428	IF ( section(is,s) >= nxl .AND. section(is,s) <= nxr ) &
1429	THEN
1430	local_2d = local_pf(section(is,s),:,nzb:nzt+1)
1431	ENDIF
1432
1433	ENDIF
1434
1435	#if defined( __parallel )
1436	IF ( netcdf_output .AND. netcdf_data_format > 4 ) THEN
1437	!
1438	!-- ATTENTION: The following lines are a workaround, because
1439	!-- independet output does not work with the
1440	!-- current netCDF4 installation. Therefore, data
1441	!-- are transferred from PEs having the cross
1442	!-- sections to other PEs along y having no cross
1443	!-- section data. Some of these data are the
1444	!-- output.
1445	!-- BEGIN WORKAROUND---------------------------------------
1446	IF ( npex /= 1 .AND. section(is,s) /= -1) THEN
1447	ALLOCATE( local_2d_l(nysg:nyng,nzb:nzt+1) )
1448	local_2d_l = 0.0
1449	IF ( section(is,s) >= nxl .AND. section(is,s) <= nxr )&
1450	THEN
1451	local_2d_l = local_2d
1452	ENDIF
1453	#if defined( __parallel )
1454	!
1455	!-- Distribute data over all PEs along x
1456	ngp = ( nyng-nysg+1 ) * ( nzt-nzb + 2 )
1457	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1458	CALL MPI_ALLREDUCE( local_2d_l(nysg,nzb), &
1459	local_2d(nysg,nzb), ngp, &
1460	MPI_REAL, MPI_SUM, comm1dx, ierr )
1461	#else
1462	local_2d = local_2d_l
1463	#endif
1464	DEALLOCATE( local_2d_l )
1465	ENDIF
1466	!-- END WORKAROUND-----------------------------------------
1467
1468	!
1469	!-- Output in netCDF4/HDF5 format.
1470	!-- Output only on those PEs where the respective cross
1471	!-- sections reside. Cross sections averaged along x are
1472	!-- output on the respective first PE along x (myidx=0).
1473	IF ( ( section(is,s) >= nxl .AND. &
1474	section(is,s) <= nxr ) .OR. &
1475	( section(is,s) == -1 .AND. myidx == 0 ) ) THEN
1476	!
1477	!-- Do not output redundant ghost point data except for the
1478	!-- boundaries of the total domain.
1479	#if defined( __netcdf )
1480	IF ( nyn == ny ) THEN
1481	nc_stat = NF90_PUT_VAR( id_set_yz(av), &
1482	id_var_do2d(av,if), &
1483	local_2d(nys:nyn+1,nzb:nzt+1), &
1484	start = (/ is, nys+1, 1, &
1485	do2d_yz_time_count(av) /), &
1486	count = (/ 1, nyn-nys+2, &
1487	nzt+2, 1 /) )
1488	ELSE
1489	nc_stat = NF90_PUT_VAR( id_set_yz(av), &
1490	id_var_do2d(av,if), &
1491	local_2d(nys:nyn,nzb:nzt+1), &
1492	start = (/ is, nys+1, 1, &
1493	do2d_yz_time_count(av) /), &
1494	count = (/ 1, nyn-nys+1, &
1495	nzt+2, 1 /) )
1496	ENDIF
1497
1498	CALL handle_netcdf_error( 'data_output_2d', 60 )
1499
1500	ELSE
1501	!
1502	!-- Output on other PEs. Only one point is output!!
1503	!-- ATTENTION: This is a workaround (see above)!!
1504	IF ( npex /= 1 ) THEN
1505	nc_stat = NF90_PUT_VAR( id_set_yz(av), &
1506	id_var_do2d(av,if), &
1507	local_2d(nys:nys,nzb:nzb), &
1508	start = (/ is, nys+1, 1, &
1509	do2d_yz_time_count(av) /), &
1510	count = (/ 1, 1, 1, 1 /) )
1511	CALL handle_netcdf_error( 'data_output_2d', 452 )
1512	ENDIF
1513	#endif
1514	ENDIF
1515
1516	ELSE
1517
1518	IF ( data_output_2d_on_each_pe ) THEN
1519	!
1520	!-- Output of partial arrays on each PE. If the cross
1521	!-- section does not reside on the PE, output special
1522	!-- index values.
1523	#if defined( __netcdf )
1524	IF ( netcdf_output .AND. myid == 0 ) THEN
1525	WRITE ( 23 ) time_since_reference_point, &
1526	do2d_yz_time_count(av), av
1527	ENDIF
1528	#endif
1529	DO i = 0, io_blocks-1
1530	IF ( i == io_group ) THEN
1531	IF ( ( section(is,s) >= nxl .AND. &
1532	section(is,s) <= nxr ) .OR. &
1533	( section(is,s) == -1 .AND. &
1534	nxl-1 == -1 ) ) &
1535	THEN
1536	WRITE (23) nysg, nyng, nzb, nzt+1
1537	WRITE (23) local_2d
1538	ELSE
1539	WRITE (23) -1, -1, -1, -1
1540	ENDIF
1541	ENDIF
1542	#if defined( __parallel )
1543	CALL MPI_BARRIER( comm2d, ierr )
1544	#endif
1545	ENDDO
1546
1547	ELSE
1548	!
1549	!-- PE0 receives partial arrays from all processors of the
1550	!-- respective cross section and outputs them. Here a
1551	!-- barrier has to be set, because otherwise
1552	!-- "-MPI- FATAL: Remote protocol queue full" may occur.
1553	CALL MPI_BARRIER( comm2d, ierr )
1554
1555	ngp = ( nyng-nysg+1 ) * ( nzt-nzb+2 )
1556	IF ( myid == 0 ) THEN
1557	!
1558	!-- Local array can be relocated directly.
1559	IF ( ( section(is,s) >= nxl .AND. &
1560	section(is,s) <= nxr ) .OR. &
1561	( section(is,s) == -1 .AND. nxl-1 == -1 ) ) &
1562	THEN
1563	total_2d(nysg:nyng,nzb:nzt+1) = local_2d
1564	ENDIF
1565	!
1566	!-- Receive data from all other PEs.
1567	DO n = 1, numprocs-1
1568	!
1569	!-- Receive index limits first, then array.
1570	!-- Index limits are received in arbitrary order from
1571	!-- the PEs.
1572	CALL MPI_RECV( ind(1), 4, MPI_INTEGER, &
1573	MPI_ANY_SOURCE, 0, comm2d, &
1574	status, ierr )
1575	!
1576	!-- Not all PEs have data for YZ-cross-section.
1577	IF ( ind(1) /= -9999 ) THEN
1578	sender = status(MPI_SOURCE)
1579	DEALLOCATE( local_2d )
1580	ALLOCATE( local_2d(ind(1):ind(2), &
1581	ind(3):ind(4)) )
1582	CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, &
1583	MPI_REAL, sender, 1, comm2d, &
1584	status, ierr )
1585	total_2d(ind(1):ind(2),ind(3):ind(4)) = &
1586	local_2d
1587	ENDIF
1588	ENDDO
1589	!
1590	!-- Output of the total cross-section.
1591	IF ( iso2d_output ) THEN
1592	WRITE (23) total_2d(0:ny+1,nzb:nzt+1)
1593	ENDIF
1594	!
1595	!-- Relocate the local array for the next loop increment
1596	DEALLOCATE( local_2d )
1597	ALLOCATE( local_2d(nysg:nyng,nzb:nzt+1) )
1598
1599	#if defined( __netcdf )
1600	IF ( netcdf_output ) THEN
1601	nc_stat = NF90_PUT_VAR( id_set_yz(av), &
1602	id_var_do2d(av,if), &
1603	total_2d(0:ny+1,nzb:nzt+1),&
1604	start = (/ is, 1, 1, do2d_yz_time_count(av) /), &
1605	count = (/ 1, ny+2, nz+2, 1 /) )
1606	CALL handle_netcdf_error( 'data_output_2d', 61 )
1607	ENDIF
1608	#endif
1609
1610	ELSE
1611	!
1612	!-- If the cross section resides on the PE, send the
1613	!-- local index limits, otherwise send -9999 to PE0.
1614	IF ( ( section(is,s) >= nxl .AND. &
1615	section(is,s) <= nxr ) .OR. &
1616	( section(is,s) == -1 .AND. nxl-1 == -1 ) ) &
1617	THEN
1618	ind(1) = nysg; ind(2) = nyng
1619	ind(3) = nzb; ind(4) = nzt+1
1620	ELSE
1621	ind(1) = -9999; ind(2) = -9999
1622	ind(3) = -9999; ind(4) = -9999
1623	ENDIF
1624	CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, &
1625	comm2d, ierr )
1626	!
1627	!-- If applicable, send data to PE0.
1628	IF ( ind(1) /= -9999 ) THEN
1629	CALL MPI_SEND( local_2d(nysg,nzb), ngp, &
1630	MPI_REAL, 0, 1, comm2d, ierr )
1631	ENDIF
1632	ENDIF
1633	!
1634	!-- A barrier has to be set, because otherwise some PEs may
1635	!-- proceed too fast so that PE0 may receive wrong data on
1636	!-- tag 0
1637	CALL MPI_BARRIER( comm2d, ierr )
1638	ENDIF
1639
1640	ENDIF
1641	#else
1642	IF ( iso2d_output ) THEN
1643	WRITE (23) local_2d(nys:nyn+1,nzb:nzt+1)
1644	ENDIF
1645	#if defined( __netcdf )
1646	IF ( netcdf_output ) THEN
1647	nc_stat = NF90_PUT_VAR( id_set_yz(av), &
1648	id_var_do2d(av,if), &
1649	local_2d(nys:nyn+1,nzb:nzt+1), &
1650	start = (/ is, 1, 1, do2d_xz_time_count(av) /), &
1651	count = (/ 1, ny+2, nz+2, 1 /) )
1652	CALL handle_netcdf_error( 'data_output_2d', 452 )
1653	ENDIF
1654	#endif
1655	#endif
1656	do2d_yz_n = do2d_yz_n + 1
1657	!
1658	!-- Write LOCAL-parameter set for ISO2D.
1659	IF ( myid == 0 .AND. iso2d_output ) THEN
1660	IF ( section(is,s) /= -1 ) THEN
1661	WRITE ( section_chr, '(''x = '',F8.2,'' m (GP '',I3, &
1662	&'')'')' &
1663	) section(is,s)*dx, section(is,s)
1664	ELSE
1665	section_chr = 'averaged along x'
1666	ENDIF
1667	IF ( av == 0 ) THEN
1668	rtext = TRIM( do2d(av,if) ) // ' t = ' // &
1669	TRIM( simulated_time_chr ) // ' ' // &
1670	TRIM( section_chr )
1671	ELSE
1672	rtext = TRIM( do2d(av,if) ) // ' averaged t = ' // &
1673	TRIM( simulated_time_chr ) // ' ' // &
1674	TRIM( section_chr )
1675	ENDIF
1676	WRITE (29,LOCAL)
1677	ENDIF
1678
1679	END SELECT
1680
1681	is = is + 1
1682	ENDDO loop1
1683
1684	ENDIF
1685
1686	if = if + 1
1687	l = MAX( 2, LEN_TRIM( do2d(av,if) ) )
1688	do2d_mode = do2d(av,if)(l-1:l)
1689
1690	ENDDO
1691
1692	!
1693	!-- Deallocate temporary arrays.
1694	IF ( ALLOCATED( level_z ) ) DEALLOCATE( level_z )
1695	DEALLOCATE( local_pf, local_2d )
1696	#if defined( __parallel )
1697	IF ( .NOT. data_output_2d_on_each_pe .AND. myid == 0 ) THEN
1698	DEALLOCATE( total_2d )
1699	ENDIF
1700	#endif
1701
1702	!
1703	!-- Close plot output file.
1704	file_id = 20 + s
1705
1706	IF ( data_output_2d_on_each_pe ) THEN
1707	DO i = 0, io_blocks-1
1708	IF ( i == io_group ) THEN
1709	CALL close_file( file_id )
1710	ENDIF
1711	#if defined( __parallel )
1712	CALL MPI_BARRIER( comm2d, ierr )
1713	#endif
1714	ENDDO
1715	ELSE
1716	IF ( myid == 0 ) CALL close_file( file_id )
1717	ENDIF
1718
1719
1720	CALL cpu_log (log_point(3),'data_output_2d','stop','nobarrier')
1721
1722	END SUBROUTINE data_output_2d

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |