Home

Context Navigation

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

Last change on this file since 1 was 1, checked in by raasch, 18 years ago
Initial repository layout and content
File size: 47.7 KB

Line
1	SUBROUTINE data_output_2d( mode, av )
2
3	!------------------------------------------------------------------------------!
4	! Actual revisions:
5	! -----------------
6	!
7	!
8	! Former revisions:
9	! -----------------
10	! $Log: data_output_2d.f90,v $
11	! Revision 1.5 2006/08/22 13:50:29 raasch
12	! xz and yz cross sections now up to nzt+1
13	!
14	! Revision 1.4 2006/06/06 07:55:18 raasch
15	! Error removed: p is assigned to to_be_resorted instead of w
16	!
17	! Revision 1.3 2006/06/02 15:18:21 raasch
18	! +argument "found" in call of routine user_data_output_2d,
19	! immediate return, if no output is requested
20	!
21	! Revision 1.2 2006/02/23 10:19:22 raasch
22	! Output of time-averaged data, output of averages along x, y, or z,
23	! output of user-defined quantities,
24	! section data are copied from local_pf to local_2d before they are output,
25	! output of particle concentration and mean radius,
26	! Former subroutine plot_2d renamed data_output_2d, pl2d.. renamed do2d..,
27	! anz renamed ngp, ebene renamed section, pl2d_.._anz renamed do2d_.._n
28	!
29	! Revision 1.1 2006/02/23 10:12:31 raasch
30	! Initial revision
31	!
32	! Revision 1.25 2005/10/20 15:41:41 raasch
33	! Error in output of netcdf yz-cross-sections removed
34	!
35	! Revision 1.24 2005/06/26 20:07:25 raasch
36	! +ql_c, ql_v, ql_vp
37	!
38	! Revision 1.23 2005/05/18 15:46:56 raasch
39	! Extensions for NetCDF output, output of vpt, qv and lwp
40	!
41	! Revision 1.22 2005/03/26 20:47:45 raasch
42	! local_2d is allocated with upper bound nzt (former nzt+1 was wrong),
43	! additional barrier in case that PE0 gathers data from the other
44	! PEs, this barrier had to be set at the end of the gathering because
45	! otherwise PE0 may receive wrong data on tag 0.
46	!
47	! Revision 1.21 2004/04/30 12:32:50 raasch
48	! Adjustments for enlarged grids (unequal subdomain sizes)
49	!
50	! Revision 1.20 2003/03/16 09:42:27 raasch
51	! Two underscores (_) are placed in front of all define-strings
52	!
53	! Revision 1.19 2002/12/19 15:56:24 raasch
54	! STOP statement replaced by call of subroutine local_stop
55	!
56	! Revision 1.18 2002/06/11 13:13:00 raasch
57	! TRIM function added (needed on IBM-systems)
58	!
59	! Revision 1.17 2002/04/16 08:11:16 08:11:16 raasch (Siegfried Raasch)
60	! Output of scalar field allowed with strings s_xy, s_xz, s_yz
61	!
62	! Revision 1.16 2001/03/30 07:44:29 raasch
63	! Translation of remaining German identifiers (variables, subroutines, etc.)
64	!
65	! Revision 1.15 2001/03/09 21:36:29 letzel
66	! All comments translated into English.
67	!
68	! Revision 1.14 2001/01/29 12:32:07 raasch
69	! Passive scalar is considered, xy-plot of q and ql corrected (level_z was
70	! missing)
71	!
72	! Revision 1.13 2001/01/22 07:38:20 raasch
73	! Module test_variables removed
74	!
75	! Revision 1.12 2001/01/02 17:33:39 raasch
76	! Plot output file is closed at the end of the subroutine
77	!
78	! Revision 1.11 2000/04/13 14:01:44 schroeter
79	! Plotausgaben fuer Gesamtwasser- und Fluessigwassergehalt
80	!
81	! Revision 1.10 99/02/05 09:14:06 09:14:06 raasch (Siegfried Raasch)
82	! u=v=0 am unteren Rand bei Vertikalschnitten wird jetzt ausschliesslich
83	! bei Spiegelungsrandbedingungen verwendet
84	!
85	! Revision 1.9 1998/09/22 17:25:58 raasch
86	! u- und theta-Ausgabe korrigiert (muessen auf nzb+1 des Hilfsfeldes local_2d
87	! umgespeichert werden.
88	!
89	! Revision 1.8 1998/07/06 12:26:06 raasch
90	! Exakte Ausgabe der Lage der Schnittebenen fuer die verschiedenen Felder
91	! bei Horizontalschnitten, + USE test_variables
92	!
93	! Revision 1.7 1998/03/09 16:24:19 raasch
94	! TRIM( simulated_time_chr )
95	!
96	! Revision 1.6 1998/02/19 07:10:30 raasch
97	! Ausgabe der Teilfelder jetzt auch wieder wahlweise durch die einzelnen PEs
98	!
99	! Revision 1.5 1998/02/10 15:08:25 raasch
100	! Spiegelungsrandbedingung fuer u und v bei Vertikalschnitten durch
101	! u=v=0 ersetzt, Ausgabe der Schnittlevel,
102	! Empfang der Teilfelder durch PE0 erfolgt in beliebiger Reihenfolge durch
103	! Verwendung von MPI_ANY_TAG, Barrieren werden gesetzt, um MPI-Fehler zu umgehen
104	!
105	! Revision 1.4 1998/01/23 10:05:42 raasch
106	! Ausgabe von 2D-Feldern (u_xy, t_xy)
107	!
108	! Revision 1.3 1997/09/08 10:47:23 raasch
109	! Prozessor 0 gibt jetzt gesamte Schnittebene des Modells aus
110	!
111	! Revision 1.2 1997/08/29 09:00:36 raasch
112	! Ausgabe der potentiellen Temperatur
113	!
114	! Revision 1.1 1997/08/11 06:24:09 raasch
115	! Initial revision
116	!
117	!
118	! Description:
119	! ------------
120	! Data output of horizontal cross-sections in NetCDF format or binary format
121	! compatible to old graphic software iso2d.
122	! Attention: The position of the sectional planes is still not always computed
123	! --------- correctly. (zu is used always)!
124	!------------------------------------------------------------------------------!
125
126	USE arrays_3d
127	USE averaging
128	USE cloud_parameters
129	USE control_parameters
130	USE cpulog
131	USE grid_variables
132	USE indices
133	USE interfaces
134	USE netcdf_control
135	USE particle_attributes
136	USE pegrid
137
138	IMPLICIT NONE
139
140	CHARACTER (LEN=2) :: do2d_mode, mode
141	CHARACTER (LEN=4) :: grid
142	CHARACTER (LEN=25) :: section_chr
143	CHARACTER (LEN=50) :: rtext
144	INTEGER :: av, ngp, file_id, i, if, is, j, k, l, layer_xy, n, psi, s, &
145	sender, &
146	ind(4)
147	LOGICAL :: found, resorted, two_d
148	REAL :: mean_r, s_r3, s_r4
149	REAL, DIMENSION(:), ALLOCATABLE :: level_z
150	REAL, DIMENSION(:,:), ALLOCATABLE :: local_2d, local_2d_l
151	REAL, DIMENSION(:,:,:), ALLOCATABLE :: local_pf
152	#if defined( __parallel )
153	REAL, DIMENSION(:,:), ALLOCATABLE :: total_2d
154	#endif
155	REAL, DIMENSION(:,:,:), POINTER :: to_be_resorted
156
157	NAMELIST /LOCAL/ rtext
158
159	CALL cpu_log (log_point(3),'data_output_2d','start')
160
161	!
162	!-- Immediate return, if no output is requested (no respective sections
163	!-- found in parameter data_output)
164	IF ( mode == 'xy' .AND. .NOT. data_output_xy(av) ) RETURN
165	IF ( mode == 'xz' .AND. .NOT. data_output_xz(av) ) RETURN
166	IF ( mode == 'yz' .AND. .NOT. data_output_yz(av) ) RETURN
167
168	two_d = .FALSE. ! local variable to distinguish between output of pure 2D
169	! arrays and cross-sections of 3D arrays.
170
171	!
172	!-- Depending on the orientation of the cross-section, the respective output
173	!-- files have to be opened.
174	SELECT CASE ( mode )
175
176	CASE ( 'xy' )
177
178	s = 1
179	ALLOCATE( level_z(0:nzt+1), local_2d(nxl-1:nxr+1,nys-1:nyn+1) )
180
181	#if defined( __netcdf )
182	IF ( myid == 0 .AND. netcdf_output ) CALL check_open( 101+av*10 )
183	#endif
184
185	IF ( data_output_2d_on_each_pe ) THEN
186	CALL check_open( 21 )
187	ELSE
188	IF ( myid == 0 ) THEN
189	IF ( iso2d_output ) CALL check_open( 21 )
190	#if defined( __parallel )
191	ALLOCATE( total_2d(-1:nx+1,-1:ny+1) )
192	#endif
193	ENDIF
194	ENDIF
195
196	CASE ( 'xz' )
197
198	s = 2
199	ALLOCATE( local_2d(nxl-1:nxr+1,nzb:nzt+1) )
200
201	#if defined( __netcdf )
202	IF ( myid == 0 .AND. netcdf_output ) CALL check_open( 102+av*10 )
203	#endif
204
205	IF ( data_output_2d_on_each_pe ) THEN
206	CALL check_open( 22 )
207	ELSE
208	IF ( myid == 0 ) THEN
209	IF ( iso2d_output ) CALL check_open( 22 )
210	#if defined( __parallel )
211	ALLOCATE( total_2d(-1:nx+1,nzb:nzt+1) )
212	#endif
213	ENDIF
214	ENDIF
215
216	CASE ( 'yz' )
217
218	s = 3
219	ALLOCATE( local_2d(nys-1:nyn+1,nzb:nzt+1) )
220
221	#if defined( __netcdf )
222	IF ( myid == 0 .AND. netcdf_output ) CALL check_open( 103+av*10 )
223	#endif
224
225	IF ( data_output_2d_on_each_pe ) THEN
226	CALL check_open( 23 )
227	ELSE
228	IF ( myid == 0 ) THEN
229	IF ( iso2d_output ) CALL check_open( 23 )
230	#if defined( __parallel )
231	ALLOCATE( total_2d(-1:ny+1,nzb:nzt+1) )
232	#endif
233	ENDIF
234	ENDIF
235
236	CASE DEFAULT
237
238	PRINT*,'+++ data_output_2d: unknown cross-section: ',mode
239	CALL local_stop
240
241	END SELECT
242
243	!
244	!-- Allocate a temporary array for resorting (kji -> ijk).
245	ALLOCATE( local_pf(nxl-1:nxr+1,nys-1:nyn+1,nzb:nzt+1) )
246
247	!
248	!-- Loop of all variables to be written.
249	!-- Output dimensions chosen
250	if = 1
251	l = MAX( 2, LEN_TRIM( do2d(av,if) ) )
252	do2d_mode = do2d(av,if)(l-1:l)
253
254	DO WHILE ( do2d(av,if)(1:1) /= ' ' )
255
256	IF ( do2d_mode == mode ) THEN
257	!
258	!-- Store the array chosen on the temporary array.
259	resorted = .FALSE.
260	SELECT CASE ( TRIM( do2d(av,if) ) )
261
262	CASE ( 'e_xy', 'e_xz', 'e_yz' )
263	IF ( av == 0 ) THEN
264	to_be_resorted => e
265	ELSE
266	to_be_resorted => e_av
267	ENDIF
268	IF ( mode == 'xy' ) level_z = zu
269
270	CASE ( 'lwp*_xy' ) ! 2d-array
271	IF ( av == 0 ) THEN
272	DO i = nxl-1, nxr+1
273	DO j = nys-1, nyn+1
274	local_pf(i,j,nzb+1) = SUM( ql(nzb:nzt,j,i) * &
275	dzw(1:nzt+1) )
276	ENDDO
277	ENDDO
278	ELSE
279	DO i = nxl-1, nxr+1
280	DO j = nys-1, nyn+1
281	local_pf(i,j,nzb+1) = lwp_av(j,i)
282	ENDDO
283	ENDDO
284	ENDIF
285	resorted = .TRUE.
286	two_d = .TRUE.
287	level_z(nzb+1) = zu(nzb+1)
288
289	CASE ( 'p_xy', 'p_xz', 'p_yz' )
290	IF ( av == 0 ) THEN
291	to_be_resorted => p
292	ELSE
293	to_be_resorted => p_av
294	ENDIF
295	IF ( mode == 'xy' ) level_z = zu
296
297	CASE ( 'pc_xy', 'pc_xz', 'pc_yz' ) ! particle concentration
298	IF ( av == 0 ) THEN
299	tend = prt_count
300	CALL exchange_horiz( tend, 0, 0 )
301	DO i = nxl-1, nxr+1
302	DO j = nys-1, nyn+1
303	DO k = nzb, nzt+1
304	local_pf(i,j,k) = tend(k,j,i)
305	ENDDO
306	ENDDO
307	ENDDO
308	resorted = .TRUE.
309	ELSE
310	CALL exchange_horiz( pc_av, 0, 0 )
311	to_be_resorted => pc_av
312	ENDIF
313
314	CASE ( 'pr_xy', 'pr_xz', 'pr_yz' ) ! mean particle radius
315	IF ( av == 0 ) THEN
316	DO i = nxl, nxr
317	DO j = nys, nyn
318	DO k = nzb, nzt+1
319	psi = prt_start_index(k,j,i)
320	s_r3 = 0.0
321	s_r4 = 0.0
322	DO n = psi, psi+prt_count(k,j,i)-1
323	s_r3 = s_r3 + particles(n)%radius**3
324	s_r4 = s_r4 + particles(n)%radius**4
325	ENDDO
326	IF ( s_r3 /= 0.0 ) THEN
327	mean_r = s_r4 / s_r3
328	ELSE
329	mean_r = 0.0
330	ENDIF
331	tend(k,j,i) = mean_r
332	ENDDO
333	ENDDO
334	ENDDO
335	CALL exchange_horiz( tend, 0, 0 )
336	DO i = nxl-1, nxr+1
337	DO j = nys-1, nyn+1
338	DO k = nzb, nzt+1
339	local_pf(i,j,k) = tend(k,j,i)
340	ENDDO
341	ENDDO
342	ENDDO
343	resorted = .TRUE.
344	ELSE
345	CALL exchange_horiz( pr_av, 0, 0 )
346	to_be_resorted => pr_av
347	ENDIF
348
349	CASE ( 'pt_xy', 'pt_xz', 'pt_yz' )
350	IF ( av == 0 ) THEN
351	IF ( .NOT. cloud_physics ) THEN
352	to_be_resorted => pt
353	ELSE
354	DO i = nxl-1, nxr+1
355	DO j = nys-1, nyn+1
356	DO k = nzb, nzt+1
357	local_pf(i,j,k) = pt(k,j,i) + l_d_cp * &
358	pt_d_t(k) * &
359	ql(k,j,i)
360	ENDDO
361	ENDDO
362	ENDDO
363	resorted = .TRUE.
364	ENDIF
365	ELSE
366	to_be_resorted => pt_av
367	ENDIF
368	IF ( mode == 'xy' ) level_z = zu
369
370	CASE ( 'q_xy', 'q_xz', 'q_yz' )
371	IF ( av == 0 ) THEN
372	to_be_resorted => q
373	ELSE
374	to_be_resorted => q_av
375	ENDIF
376	IF ( mode == 'xy' ) level_z = zu
377
378	CASE ( 'ql_xy', 'ql_xz', 'ql_yz' )
379	IF ( av == 0 ) THEN
380	to_be_resorted => ql
381	ELSE
382	to_be_resorted => ql_av
383	ENDIF
384	IF ( mode == 'xy' ) level_z = zu
385
386	CASE ( 'ql_c_xy', 'ql_c_xz', 'ql_c_yz' )
387	IF ( av == 0 ) THEN
388	to_be_resorted => ql_c
389	ELSE
390	to_be_resorted => ql_c_av
391	ENDIF
392	IF ( mode == 'xy' ) level_z = zu
393
394	CASE ( 'ql_v_xy', 'ql_v_xz', 'ql_v_yz' )
395	IF ( av == 0 ) THEN
396	to_be_resorted => ql_v
397	ELSE
398	to_be_resorted => ql_v_av
399	ENDIF
400	IF ( mode == 'xy' ) level_z = zu
401
402	CASE ( 'ql_vp_xy', 'ql_vp_xz', 'ql_vp_yz' )
403	IF ( av == 0 ) THEN
404	to_be_resorted => ql_vp
405	ELSE
406	to_be_resorted => ql_vp_av
407	ENDIF
408	IF ( mode == 'xy' ) level_z = zu
409
410	CASE ( 'qv_xy', 'qv_xz', 'qv_yz' )
411	IF ( av == 0 ) THEN
412	DO i = nxl-1, nxr+1
413	DO j = nys-1, nyn+1
414	DO k = nzb, nzt+1
415	local_pf(i,j,k) = q(k,j,i) - ql(k,j,i)
416	ENDDO
417	ENDDO
418	ENDDO
419	resorted = .TRUE.
420	ELSE
421	to_be_resorted => qv_av
422	ENDIF
423	IF ( mode == 'xy' ) level_z = zu
424
425	CASE ( 's_xy', 's_xz', 's_yz' )
426	IF ( av == 0 ) THEN
427	to_be_resorted => q
428	ELSE
429	to_be_resorted => q_av
430	ENDIF
431
432	CASE ( 't*_xy' ) ! 2d-array
433	IF ( av == 0 ) THEN
434	DO i = nxl-1, nxr+1
435	DO j = nys-1, nyn+1
436	local_pf(i,j,nzb+1) = ts(j,i)
437	ENDDO
438	ENDDO
439	ELSE
440	DO i = nxl-1, nxr+1
441	DO j = nys-1, nyn+1
442	local_pf(i,j,nzb+1) = ts_av(j,i)
443	ENDDO
444	ENDDO
445	ENDIF
446	resorted = .TRUE.
447	two_d = .TRUE.
448	level_z(nzb+1) = zu(nzb+1)
449
450	CASE ( 'u_xy', 'u_xz', 'u_yz' )
451	IF ( av == 0 ) THEN
452	to_be_resorted => u
453	ELSE
454	to_be_resorted => u_av
455	ENDIF
456	IF ( mode == 'xy' ) level_z = zu
457	!
458	!-- Substitute the values generated by "mirror" boundary condition
459	!-- at the bottom boundary by the real surface values.
460	IF ( do2d(av,if) == 'u_xz' .OR. do2d(av,if) == 'u_yz' ) THEN
461	IF ( ibc_uv_b == 0 ) local_pf(:,:,nzb) = 0.0
462	ENDIF
463
464	CASE ( 'u*_xy' ) ! 2d-array
465	IF ( av == 0 ) THEN
466	DO i = nxl-1, nxr+1
467	DO j = nys-1, nyn+1
468	local_pf(i,j,nzb+1) = us(j,i)
469	ENDDO
470	ENDDO
471	ELSE
472	DO i = nxl-1, nxr+1
473	DO j = nys-1, nyn+1
474	local_pf(i,j,nzb+1) = us_av(j,i)
475	ENDDO
476	ENDDO
477	ENDIF
478	resorted = .TRUE.
479	two_d = .TRUE.
480	level_z(nzb+1) = zu(nzb+1)
481
482	CASE ( 'v_xy', 'v_xz', 'v_yz' )
483	IF ( av == 0 ) THEN
484	to_be_resorted => v
485	ELSE
486	to_be_resorted => v_av
487	ENDIF
488	IF ( mode == 'xy' ) level_z = zu
489	!
490	!-- Substitute the values generated by "mirror" boundary condition
491	!-- at the bottom boundary by the real surface values.
492	IF ( do2d(av,if) == 'v_xz' .OR. do2d(av,if) == 'v_yz' ) THEN
493	IF ( ibc_uv_b == 0 ) local_pf(:,:,nzb) = 0.0
494	ENDIF
495
496	CASE ( 'vpt_xy', 'vpt_xz', 'vpt_yz' )
497	IF ( av == 0 ) THEN
498	to_be_resorted => vpt
499	ELSE
500	to_be_resorted => vpt_av
501	ENDIF
502	IF ( mode == 'xy' ) level_z = zu
503
504	CASE ( 'w_xy', 'w_xz', 'w_yz' )
505	IF ( av == 0 ) THEN
506	to_be_resorted => w
507	ELSE
508	to_be_resorted => w_av
509	ENDIF
510	IF ( mode == 'xy' ) level_z = zw
511
512	CASE DEFAULT
513	!
514	!-- User defined quantity
515	CALL user_data_output_2d( av, do2d(av,if), found, grid, &
516	local_pf )
517	resorted = .TRUE.
518
519	IF ( grid == 'zu' ) THEN
520	IF ( mode == 'xy' ) level_z = zu
521	ELSEIF ( grid == 'zw' ) THEN
522	IF ( mode == 'xy' ) level_z = zw
523	ENDIF
524
525	IF ( .NOT. found ) THEN
526	PRINT*, '+++ data_output_2d: no output provided for: ', &
527	do2d(av,if)
528	ENDIF
529
530	END SELECT
531
532	!
533	!-- Resort the array to be output, if not done above
534	IF ( .NOT. resorted ) THEN
535	DO i = nxl-1, nxr+1
536	DO j = nys-1, nyn+1
537	DO k = nzb, nzt+1
538	local_pf(i,j,k) = to_be_resorted(k,j,i)
539	ENDDO
540	ENDDO
541	ENDDO
542	ENDIF
543
544	!
545	!-- Output of the individual cross-sections, depending on the cross-
546	!-- section mode chosen.
547	is = 1
548	loop1: DO WHILE ( section(is,s) /= -9999 .OR. two_d )
549
550	SELECT CASE ( mode )
551
552	CASE ( 'xy' )
553	!
554	!-- Determine the cross section index
555	IF ( two_d ) THEN
556	layer_xy = nzb+1
557	ELSE
558	layer_xy = section(is,s)
559	ENDIF
560
561	!
562	!-- Update the NetCDF xy cross section time axis
563	IF ( myid == 0 ) THEN
564	IF ( simulated_time /= do2d_xy_last_time(av) ) THEN
565	do2d_xy_time_count(av) = do2d_xy_time_count(av) + 1
566	do2d_xy_last_time(av) = simulated_time
567	IF ( .NOT. data_output_2d_on_each_pe .AND. &
568	netcdf_output ) THEN
569	#if defined( __netcdf )
570	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
571	id_var_time_xy(av), &
572	(/ simulated_time /), &
573	start = (/ do2d_xy_time_count(av) /), &
574	count = (/ 1 /) )
575	IF ( nc_stat /= NF90_NOERR ) THEN
576	CALL handle_netcdf_error( 53 )
577	ENDIF
578	#endif
579	ENDIF
580	ENDIF
581	ENDIF
582	!
583	!-- If required, carry out averaging along z
584	IF ( section(is,s) == -1 ) THEN
585
586	local_2d = 0.0
587	!
588	!-- Carry out the averaging (all data are on the PE)
589	DO k = nzb, nzt+1
590	DO j = nys-1, nyn+1
591	DO i = nxl-1, nxr+1
592	local_2d(i,j) = local_2d(i,j) + local_pf(i,j,k)
593	ENDDO
594	ENDDO
595	ENDDO
596
597	local_2d = local_2d / ( nzt -nzb + 2.0 )
598
599	ELSE
600	!
601	!-- Just store the respective section on the local array
602	local_2d = local_pf(:,:,layer_xy)
603
604	ENDIF
605
606	#if defined( __parallel )
607	IF ( data_output_2d_on_each_pe ) THEN
608	!
609	!-- Output of partial arrays on each PE
610	#if defined( __netcdf )
611	IF ( netcdf_output .AND. myid == 0 ) THEN
612	WRITE ( 21 ) simulated_time, do2d_xy_time_count(av), &
613	av
614	ENDIF
615	#endif
616	WRITE ( 21 ) nxl-1, nxr+1, nys-1, nyn+1
617	WRITE ( 21 ) local_2d
618
619	ELSE
620	!
621	!-- PE0 receives partial arrays from all processors and then
622	!-- outputs them. Here a barrier has to be set, because
623	!-- otherwise "-MPI- FATAL: Remote protocol queue full" may
624	!-- occur.
625	CALL MPI_BARRIER( comm2d, ierr )
626
627	ngp = ( nxr-nxl+3 ) * ( nyn-nys+3 )
628	IF ( myid == 0 ) THEN
629	!
630	!-- Local array can be relocated directly.
631	total_2d(nxl-1:nxr+1,nys-1:nyn+1) = local_2d
632	!
633	!-- Receive data from all other PEs.
634	DO n = 1, numprocs-1
635	!
636	!-- Receive index limits first, then array.
637	!-- Index limits are received in arbitrary order from
638	!-- the PEs.
639	CALL MPI_RECV( ind(1), 4, MPI_INTEGER, &
640	MPI_ANY_SOURCE, 0, comm2d, status, &
641	ierr )
642	sender = status(MPI_SOURCE)
643	DEALLOCATE( local_2d )
644	ALLOCATE( local_2d(ind(1):ind(2),ind(3):ind(4)) )
645	CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, &
646	MPI_REAL, sender, 1, comm2d, &
647	status, ierr )
648	total_2d(ind(1):ind(2),ind(3):ind(4)) = local_2d
649	ENDDO
650	!
651	!-- Output of the total cross-section.
652	IF ( iso2d_output ) WRITE (21) total_2d(0:nx+1,0:ny+1)
653	!
654	!-- Relocate the local array for the next loop increment
655	DEALLOCATE( local_2d )
656	ALLOCATE( local_2d(nxl-1:nxr+1,nys-1:nyn+1) )
657
658	#if defined( __netcdf )
659	IF ( netcdf_output ) THEN
660	IF ( two_d ) THEN
661	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
662	id_var_do2d(av,if), &
663	total_2d(0:nx+1,0:ny+1), &
664	start = (/ 1, 1, 1, do2d_xy_time_count(av) /), &
665	count = (/ nx+2, ny+2, 1, 1 /) )
666	ELSE
667	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
668	id_var_do2d(av,if), &
669	total_2d(0:nx+1,0:ny+1), &
670	start = (/ 1, 1, is, do2d_xy_time_count(av) /), &
671	count = (/ nx+2, ny+2, 1, 1 /) )
672	ENDIF
673	IF ( nc_stat /= NF90_NOERR ) &
674	CALL handle_netcdf_error( 54 )
675	ENDIF
676	#endif
677
678	ELSE
679	!
680	!-- First send the local index limits to PE0
681	ind(1) = nxl-1; ind(2) = nxr+1
682	ind(3) = nys-1; ind(4) = nyn+1
683	CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, comm2d, &
684	ierr )
685	!
686	!-- Send data to PE0
687	CALL MPI_SEND( local_2d(nxl-1,nys-1), ngp, MPI_REAL, &
688	0, 1, comm2d, ierr )
689	ENDIF
690	!
691	!-- A barrier has to be set, because otherwise some PEs may
692	!-- proceed too fast so that PE0 may receive wrong data on
693	!-- tag 0
694	CALL MPI_BARRIER( comm2d, ierr )
695	ENDIF
696	#else
697	IF ( iso2d_output ) THEN
698	WRITE (21) local_2d(nxl:nxr+1,nys:nyn+1)
699	ENDIF
700	#if defined( __netcdf )
701	IF ( netcdf_output ) THEN
702	IF ( two_d ) THEN
703	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
704	id_var_do2d(av,if), &
705	local_2d(nxl:nxr+1,nys:nyn+1), &
706	start = (/ 1, 1, 1, do2d_xy_time_count(av) /), &
707	count = (/ nx+2, ny+2, 1, 1 /) )
708	ELSE
709	nc_stat = NF90_PUT_VAR( id_set_xy(av), &
710	id_var_do2d(av,if), &
711	local_2d(nxl:nxr+1,nys:nyn+1), &
712	start = (/ 1, 1, is, do2d_xy_time_count(av) /), &
713	count = (/ nx+2, ny+2, 1, 1 /) )
714	ENDIF
715	IF ( nc_stat /= NF90_NOERR ) &
716	CALL handle_netcdf_error( 55 )
717	ENDIF
718	#endif
719	#endif
720	do2d_xy_n = do2d_xy_n + 1
721	!
722	!-- Write LOCAL parameter set for ISO2D.
723	IF ( myid == 0 .AND. iso2d_output ) THEN
724	IF ( section(is,s) /= -1 ) THEN
725	WRITE ( section_chr, '(''z = '',F7.2,'' m (GP '',I3, &
726	&'')'')' &
727	) level_z(layer_xy), layer_xy
728	ELSE
729	section_chr = 'averaged along z'
730	ENDIF
731	IF ( av == 0 ) THEN
732	rtext = TRIM( do2d(av,if) ) // ' t = ' // &
733	TRIM( simulated_time_chr ) // ' ' // &
734	TRIM( section_chr )
735	ELSE
736	rtext = TRIM( do2d(av,if) ) // ' averaged t = ' // &
737	TRIM( simulated_time_chr ) // ' ' // &
738	TRIM( section_chr )
739	ENDIF
740	WRITE (27,LOCAL)
741	ENDIF
742	!
743	!-- For 2D-arrays (e.g. u*) only one cross-section is available.
744	!-- Hence exit loop of output levels.
745	IF ( two_d ) THEN
746	two_d = .FALSE.
747	EXIT loop1
748	ENDIF
749
750	CASE ( 'xz' )
751	!
752	!-- Update the NetCDF xy cross section time axis
753	IF ( myid == 0 ) THEN
754	IF ( simulated_time /= do2d_xz_last_time(av) ) THEN
755	do2d_xz_time_count(av) = do2d_xz_time_count(av) + 1
756	do2d_xz_last_time(av) = simulated_time
757	IF ( .NOT. data_output_2d_on_each_pe .AND. &
758	netcdf_output ) THEN
759	#if defined( __netcdf )
760	nc_stat = NF90_PUT_VAR( id_set_xz(av), &
761	id_var_time_xz(av), &
762	(/ simulated_time /), &
763	start = (/ do2d_xz_time_count(av) /), &
764	count = (/ 1 /) )
765	IF ( nc_stat /= NF90_NOERR ) THEN
766	CALL handle_netcdf_error( 56 )
767	ENDIF
768	#endif
769	ENDIF
770	ENDIF
771	ENDIF
772	!
773	!-- If required, carry out averaging along y
774	IF ( section(is,s) == -1 ) THEN
775
776	ALLOCATE( local_2d_l(nxl-1:nxr+1,nzb:nzt+1) )
777	local_2d_l = 0.0
778	ngp = ( nxr-nxl+3 ) * ( nzt-nzb+2 )
779	!
780	!-- First local averaging on the PE
781	DO k = nzb, nzt+1
782	DO j = nys, nyn
783	DO i = nxl-1, nxr+1
784	local_2d_l(i,k) = local_2d_l(i,k) + &
785	local_pf(i,j,k)
786	ENDDO
787	ENDDO
788	ENDDO
789	#if defined( __parallel )
790	!
791	!-- Now do the averaging over all PEs along y
792	CALL MPI_ALLREDUCE( local_2d_l(nxl-1,nzb), &
793	local_2d(nxl-1,nzb), ngp, MPI_REAL, &
794	MPI_SUM, comm1dy, ierr )
795	#else
796	local_2d = local_2d_l
797	#endif
798	local_2d = local_2d / ( ny + 1.0 )
799
800	DEALLOCATE( local_2d_l )
801
802	ELSE
803	!
804	!-- Just store the respective section on the local array
805	!-- (but only if it is available on this PE!)
806	IF ( section(is,s) >= nys .AND. section(is,s) <= nyn ) &
807	THEN
808	local_2d = local_pf(:,section(is,s),nzb:nzt+1)
809	ENDIF
810
811	ENDIF
812
813	#if defined( __parallel )
814	IF ( data_output_2d_on_each_pe ) THEN
815	!
816	!-- Output of partial arrays on each PE. If the cross section
817	!-- does not reside on the PE, output special index values.
818	#if defined( __netcdf )
819	IF ( netcdf_output .AND. myid == 0 ) THEN
820	WRITE ( 22 ) simulated_time, do2d_xz_time_count(av), &
821	av
822	ENDIF
823	#endif
824	IF ( ( section(is,s)>=nys .AND. section(is,s)<=nyn ) .OR.&
825	( section(is,s) == -1 .AND. nys-1 == -1 ) ) &
826	THEN
827	WRITE (22) nxl-1, nxr+1, nzb, nzt+1
828	WRITE (22) local_2d
829	ELSE
830	WRITE (22) -1, -1, -1, -1
831	ENDIF
832
833	ELSE
834	!
835	!-- PE0 receives partial arrays from all processors of the
836	!-- respective cross section and outputs them. Here a
837	!-- barrier has to be set, because otherwise
838	!-- "-MPI- FATAL: Remote protocol queue full" may occur.
839	CALL MPI_BARRIER( comm2d, ierr )
840
841	ngp = ( nxr-nxl+3 ) * ( nzt-nzb+2 )
842	IF ( myid == 0 ) THEN
843	!
844	!-- Local array can be relocated directly.
845	IF ( ( section(is,s)>=nys .AND. section(is,s)<=nyn ) &
846	.OR. ( section(is,s) == -1 .AND. nys-1 == -1 ) ) &
847	THEN
848	total_2d(nxl-1:nxr+1,nzb:nzt+1) = local_2d
849	ENDIF
850	!
851	!-- Receive data from all other PEs.
852	DO n = 1, numprocs-1
853	!
854	!-- Receive index limits first, then array.
855	!-- Index limits are received in arbitrary order from
856	!-- the PEs.
857	CALL MPI_RECV( ind(1), 4, MPI_INTEGER, &
858	MPI_ANY_SOURCE, 0, comm2d, status, &
859	ierr )
860	!
861	!-- Not all PEs have data for XZ-cross-section.
862	IF ( ind(1) /= -9999 ) THEN
863	sender = status(MPI_SOURCE)
864	DEALLOCATE( local_2d )
865	ALLOCATE( local_2d(ind(1):ind(2),ind(3):ind(4)) )
866	CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, &
867	MPI_REAL, sender, 1, comm2d, &
868	status, ierr )
869	total_2d(ind(1):ind(2),ind(3):ind(4)) = local_2d
870	ENDIF
871	ENDDO
872	!
873	!-- Output of the total cross-section.
874	IF ( iso2d_output ) THEN
875	WRITE (22) total_2d(0:nx+1,nzb:nzt+1)
876	ENDIF
877	!
878	!-- Relocate the local array for the next loop increment
879	DEALLOCATE( local_2d )
880	ALLOCATE( local_2d(nxl-1:nxr+1,nzb:nzt+1) )
881
882	#if defined( __netcdf )
883	IF ( netcdf_output ) THEN
884	nc_stat = NF90_PUT_VAR( id_set_xz(av), &
885	id_var_do2d(av,if), &
886	total_2d(0:nx+1,nzb:nzt+1),&
887	start = (/ 1, is, 1, do2d_xz_time_count(av) /), &
888	count = (/ nx+2, 1, nz+2, 1 /) )
889	IF ( nc_stat /= NF90_NOERR ) &
890	CALL handle_netcdf_error( 57 )
891	ENDIF
892	#endif
893
894	ELSE
895	!
896	!-- If the cross section resides on the PE, send the
897	!-- local index limits, otherwise send -9999 to PE0.
898	IF ( ( section(is,s)>=nys .AND. section(is,s)<=nyn ) &
899	.OR. ( section(is,s) == -1 .AND. nys-1 == -1 ) ) &
900	THEN
901	ind(1) = nxl-1; ind(2) = nxr+1
902	ind(3) = nzb; ind(4) = nzt+1
903	ELSE
904	ind(1) = -9999; ind(2) = -9999
905	ind(3) = -9999; ind(4) = -9999
906	ENDIF
907	CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, comm2d, &
908	ierr )
909	!
910	!-- If applicable, send data to PE0.
911	IF ( ind(1) /= -9999 ) THEN
912	CALL MPI_SEND( local_2d(nxl-1,nzb), ngp, MPI_REAL, &
913	0, 1, comm2d, ierr )
914	ENDIF
915	ENDIF
916	!
917	!-- A barrier has to be set, because otherwise some PEs may
918	!-- proceed too fast so that PE0 may receive wrong data on
919	!-- tag 0
920	CALL MPI_BARRIER( comm2d, ierr )
921	ENDIF
922	#else
923	IF ( iso2d_output ) THEN
924	WRITE (22) local_2d(nxl:nxr+1,nzb:nzt+1)
925	ENDIF
926	#if defined( __netcdf )
927	IF ( netcdf_output ) THEN
928	nc_stat = NF90_PUT_VAR( id_set_xz(av), &
929	id_var_do2d(av,if), &
930	local_2d(nxl:nxr+1,nzb:nzt+1), &
931	start = (/ 1, is, 1, do2d_xz_time_count(av) /), &
932	count = (/ nx+2, 1, nz+2, 1 /) )
933	IF ( nc_stat /= NF90_NOERR ) &
934	CALL handle_netcdf_error( 58 )
935	ENDIF
936	#endif
937	#endif
938	do2d_xz_n = do2d_xz_n + 1
939	!
940	!-- Write LOCAL-parameter set for ISO2D.
941	IF ( myid == 0 .AND. iso2d_output ) THEN
942	IF ( section(is,s) /= -1 ) THEN
943	WRITE ( section_chr, '(''y = '',F8.2,'' m (GP '',I3, &
944	&'')'')' &
945	) section(is,s)*dy, section(is,s)
946	ELSE
947	section_chr = 'averaged along y'
948	ENDIF
949	IF ( av == 0 ) THEN
950	rtext = TRIM( do2d(av,if) ) // ' t = ' // &
951	TRIM( simulated_time_chr ) // ' ' // &
952	TRIM( section_chr )
953	ELSE
954	rtext = TRIM( do2d(av,if) ) // ' averaged t = ' // &
955	TRIM( simulated_time_chr ) // ' ' // &
956	TRIM( section_chr )
957	ENDIF
958	WRITE (28,LOCAL)
959	ENDIF
960
961	CASE ( 'yz' )
962	!
963	!-- Update the NetCDF xy cross section time axis
964	IF ( myid == 0 ) THEN
965	IF ( simulated_time /= do2d_yz_last_time(av) ) THEN
966	do2d_yz_time_count(av) = do2d_yz_time_count(av) + 1
967	do2d_yz_last_time(av) = simulated_time
968	IF ( .NOT. data_output_2d_on_each_pe .AND. &
969	netcdf_output ) THEN
970	#if defined( __netcdf )
971	nc_stat = NF90_PUT_VAR( id_set_yz(av), &
972	id_var_time_yz(av), &
973	(/ simulated_time /), &
974	start = (/ do2d_yz_time_count(av) /), &
975	count = (/ 1 /) )
976	IF ( nc_stat /= NF90_NOERR ) THEN
977	CALL handle_netcdf_error( 59 )
978	ENDIF
979	#endif
980	ENDIF
981	ENDIF
982	ENDIF
983	!
984	!-- If required, carry out averaging along x
985	IF ( section(is,s) == -1 ) THEN
986
987	ALLOCATE( local_2d_l(nys-1:nyn+1,nzb:nzt+1) )
988	local_2d_l = 0.0
989	ngp = ( nyn-nys+3 ) * ( nzt-nzb+2 )
990	!
991	!-- First local averaging on the PE
992	DO k = nzb, nzt+1
993	DO j = nys-1, nyn+1
994	DO i = nxl, nxr
995	local_2d_l(j,k) = local_2d_l(j,k) + &
996	local_pf(i,j,k)
997	ENDDO
998	ENDDO
999	ENDDO
1000	#if defined( __parallel )
1001	!
1002	!-- Now do the averaging over all PEs along x
1003	CALL MPI_ALLREDUCE( local_2d_l(nys-1,nzb), &
1004	local_2d(nys-1,nzb), ngp, MPI_REAL, &
1005	MPI_SUM, comm1dx, ierr )
1006	#else
1007	local_2d = local_2d_l
1008	#endif
1009	local_2d = local_2d / ( nx + 1.0 )
1010
1011	DEALLOCATE( local_2d_l )
1012
1013	ELSE
1014	!
1015	!-- Just store the respective section on the local array
1016	!-- (but only if it is available on this PE!)
1017	IF ( section(is,s) >= nxl .AND. section(is,s) <= nxr ) &
1018	THEN
1019	local_2d = local_pf(section(is,s),:,nzb:nzt+1)
1020	ENDIF
1021
1022	ENDIF
1023
1024	#if defined( __parallel )
1025	IF ( data_output_2d_on_each_pe ) THEN
1026	!
1027	!-- Output of partial arrays on each PE. If the cross section
1028	!-- does not reside on the PE, output special index values.
1029	#if defined( __netcdf )
1030	IF ( netcdf_output .AND. myid == 0 ) THEN
1031	WRITE ( 23 ) simulated_time, do2d_yz_time_count(av), &
1032	av
1033	ENDIF
1034	#endif
1035	IF ( ( section(is,s)>=nxl .AND. section(is,s)<=nxr ) .OR.&
1036	( section(is,s) == -1 .AND. nxl-1 == -1 ) ) &
1037	THEN
1038	WRITE (23) nys-1, nyn+1, nzb, nzt+1
1039	WRITE (23) local_2d
1040	ELSE
1041	WRITE (23) -1, -1, -1, -1
1042	ENDIF
1043
1044	ELSE
1045	!
1046	!-- PE0 receives partial arrays from all processors of the
1047	!-- respective cross section and outputs them. Here a
1048	!-- barrier has to be set, because otherwise
1049	!-- "-MPI- FATAL: Remote protocol queue full" may occur.
1050	CALL MPI_BARRIER( comm2d, ierr )
1051
1052	ngp = ( nyn-nys+3 ) * ( nzt-nzb+2 )
1053	IF ( myid == 0 ) THEN
1054	!
1055	!-- Local array can be relocated directly.
1056	IF ( ( section(is,s)>=nxl .AND. section(is,s)<=nxr ) &
1057	.OR. ( section(is,s) == -1 .AND. nxl-1 == -1 ) ) &
1058	THEN
1059	total_2d(nys-1:nyn+1,nzb:nzt+1) = local_2d
1060	ENDIF
1061	!
1062	!-- Receive data from all other PEs.
1063	DO n = 1, numprocs-1
1064	!
1065	!-- Receive index limits first, then array.
1066	!-- Index limits are received in arbitrary order from
1067	!-- the PEs.
1068	CALL MPI_RECV( ind(1), 4, MPI_INTEGER, &
1069	MPI_ANY_SOURCE, 0, comm2d, status, &
1070	ierr )
1071	!
1072	!-- Not all PEs have data for YZ-cross-section.
1073	IF ( ind(1) /= -9999 ) THEN
1074	sender = status(MPI_SOURCE)
1075	DEALLOCATE( local_2d )
1076	ALLOCATE( local_2d(ind(1):ind(2),ind(3):ind(4)) )
1077	CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, &
1078	MPI_REAL, sender, 1, comm2d, &
1079	status, ierr )
1080	total_2d(ind(1):ind(2),ind(3):ind(4)) = local_2d
1081	ENDIF
1082	ENDDO
1083	!
1084	!-- Output of the total cross-section.
1085	IF ( iso2d_output ) THEN
1086	WRITE (23) total_2d(0:ny+1,nzb:nzt+1)
1087	ENDIF
1088	!
1089	!-- Relocate the local array for the next loop increment
1090	DEALLOCATE( local_2d )
1091	ALLOCATE( local_2d(nys-1:nyn+1,nzb:nzt+1) )
1092
1093	#if defined( __netcdf )
1094	IF ( netcdf_output ) THEN
1095	nc_stat = NF90_PUT_VAR( id_set_yz(av), &
1096	id_var_do2d(av,if), &
1097	total_2d(0:ny+1,nzb:nzt+1),&
1098	start = (/ is, 1, 1, do2d_yz_time_count(av) /), &
1099	count = (/ 1, ny+2, nz+2, 1 /) )
1100	IF ( nc_stat /= NF90_NOERR ) &
1101	CALL handle_netcdf_error( 60 )
1102	ENDIF
1103	#endif
1104
1105	ELSE
1106	!
1107	!-- If the cross section resides on the PE, send the
1108	!-- local index limits, otherwise send -9999 to PE0.
1109	IF ( ( section(is,s)>=nxl .AND. section(is,s)<=nxr ) &
1110	.OR. ( section(is,s) == -1 .AND. nxl-1 == -1 ) ) &
1111	THEN
1112	ind(1) = nys-1; ind(2) = nyn+1
1113	ind(3) = nzb; ind(4) = nzt+1
1114	ELSE
1115	ind(1) = -9999; ind(2) = -9999
1116	ind(3) = -9999; ind(4) = -9999
1117	ENDIF
1118	CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, comm2d, &
1119	ierr )
1120	!
1121	!-- If applicable, send data to PE0.
1122	IF ( ind(1) /= -9999 ) THEN
1123	CALL MPI_SEND( local_2d(nys-1,nzb), ngp, MPI_REAL, &
1124	0, 1, comm2d, ierr )
1125	ENDIF
1126	ENDIF
1127	!
1128	!-- A barrier has to be set, because otherwise some PEs may
1129	!-- proceed too fast so that PE0 may receive wrong data on
1130	!-- tag 0
1131	CALL MPI_BARRIER( comm2d, ierr )
1132	ENDIF
1133	#else
1134	IF ( iso2d_output ) THEN
1135	WRITE (23) local_2d(nys:nyn+1,nzb:nzt+1)
1136	ENDIF
1137	#if defined( __netcdf )
1138	IF ( netcdf_output ) THEN
1139	nc_stat = NF90_PUT_VAR( id_set_yz(av), &
1140	id_var_do2d(av,if), &
1141	local_2d(nys:nyn+1,nzb:nzt+1), &
1142	start = (/ is, 1, 1, do2d_xz_time_count(av) /), &
1143	count = (/ 1, ny+2, nz+2, 1 /) )
1144	IF ( nc_stat /= NF90_NOERR ) &
1145	CALL handle_netcdf_error( 61 )
1146	ENDIF
1147	#endif
1148	#endif
1149	do2d_yz_n = do2d_yz_n + 1
1150	!
1151	!-- Write LOCAL-parameter set for ISO2D.
1152	IF ( myid == 0 .AND. iso2d_output ) THEN
1153	IF ( section(is,s) /= -1 ) THEN
1154	WRITE ( section_chr, '(''x = '',F8.2,'' m (GP '',I3, &
1155	&'')'')' &
1156	) section(is,s)*dx, section(is,s)
1157	ELSE
1158	section_chr = 'averaged along x'
1159	ENDIF
1160	IF ( av == 0 ) THEN
1161	rtext = TRIM( do2d(av,if) ) // ' t = ' // &
1162	TRIM( simulated_time_chr ) // ' ' // &
1163	TRIM( section_chr )
1164	ELSE
1165	rtext = TRIM( do2d(av,if) ) // ' averaged t = ' // &
1166	TRIM( simulated_time_chr ) // ' ' // &
1167	TRIM( section_chr )
1168	ENDIF
1169	WRITE (29,LOCAL)
1170	ENDIF
1171
1172	END SELECT
1173
1174	is = is + 1
1175	ENDDO loop1
1176
1177	ENDIF
1178
1179	if = if + 1
1180	l = MAX( 2, LEN_TRIM( do2d(av,if) ) )
1181	do2d_mode = do2d(av,if)(l-1:l)
1182
1183	ENDDO
1184
1185	!
1186	!-- Deallocate temporary arrays.
1187	IF ( ALLOCATED( level_z ) ) DEALLOCATE( level_z )
1188	DEALLOCATE( local_pf, local_2d )
1189	#if defined( __parallel )
1190	IF ( .NOT. data_output_2d_on_each_pe .AND. myid == 0 ) THEN
1191	DEALLOCATE( total_2d )
1192	ENDIF
1193	#endif
1194
1195	!
1196	!-- Close plot output file.
1197	file_id = 20 + s
1198
1199	IF ( data_output_2d_on_each_pe ) THEN
1200	CALL close_file( file_id )
1201	ELSE
1202	IF ( myid == 0 ) CALL close_file( file_id )
1203	ENDIF
1204
1205
1206	CALL cpu_log (log_point(3),'data_output_2d','stop','nobarrier')
1207
1208	END SUBROUTINE data_output_2d

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |