Home

Context Navigation

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

Last change on this file since 181 was 108, checked in by letzel, 17 years ago

Improved coupler: evaporation - salinity-flux coupling for humidity = .T.,

avoid MPI hangs when coupled runs terminate, add DOC/app/chapter_3.8;

Optional calculation of km and kh from initial TKE e_init;
Default initialization of km,kh = 0.00001 for ocean = .T.;
Allow data_output_pr= q, wq, w"q", w*q* for humidity = .T.;
Bugfix: Rayleigh damping for ocean fixed.

Property svn:keywords set to Id

File size: 46.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |