Home

Context Navigation

source: palm/trunk/SOURCE/data_output_3d.f90 @ 1308

Last change on this file since 1308 was 1308, checked in by fricke, 10 years ago
Adjustments for parallel NetCDF output for lccrayh/lccrayb (Cray XC30 systems)
Property svn:keywords set to `Id`
File size: 21.1 KB

Line
1	SUBROUTINE data_output_3d( 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	! Check, if the limit of the time dimension is exceeded for parallel output
23	! To increase the performance for parallel output, the following is done:
24	! - Update of time axis is only done by PE0
25	!
26	! Former revisions:
27	! -----------------
28	! $Id: data_output_3d.f90 1308 2014-03-13 14:58:42Z fricke $
29	!
30	! 1244 2013-10-31 08:16:56Z raasch
31	! Bugfix for index bounds in case of 3d-parallel output
32	!
33	! 1115 2013-03-26 18:16:16Z hoffmann
34	! ql is calculated by calc_liquid_water_content
35	!
36	! 1106 2013-03-04 05:31:38Z raasch
37	! array_kind renamed precision_kind
38	!
39	! 1076 2012-12-05 08:30:18Z hoffmann
40	! Bugfix in output of ql
41	!
42	! 1053 2012-11-13 17:11:03Z hoffmann
43	! +nr, qr, prr, qc and averaged quantities
44	!
45	! 1036 2012-10-22 13:43:42Z raasch
46	! code put under GPL (PALM 3.9)
47	!
48	! 1031 2012-10-19 14:35:30Z raasch
49	! netCDF4 without parallel file support implemented
50	!
51	! 1007 2012-09-19 14:30:36Z franke
52	! Bugfix: missing calculation of ql_vp added
53	!
54	! 790 2011-11-29 03:11:20Z raasch
55	! bugfix: calculation of 'pr' must depend on the particle weighting factor,
56	! nzt+1 replaced by nz_do3d for 'pr'
57	!
58	! 771 2011-10-27 10:56:21Z heinze
59	! +lpt
60	!
61	! 759 2011-09-15 13:58:31Z raasch
62	! Splitting of parallel I/O
63	!
64	! 727 2011-04-20 20:05:25Z suehring
65	! Exchange ghost layers also for p_av.
66	!
67	! 725 2011-04-11 09:37:01Z suehring
68	! Exchange ghost layers for p regardless of used pressure solver (except SOR).
69	!
70	! 691 2011-03-04 08:45:30Z maronga
71	! Replaced simulated_time by time_since_reference_point
72	!
73	! 673 2011-01-18 16:19:48Z suehring
74	! When using Multigrid or SOR solver an additional CALL exchange_horiz is
75	! is needed for pressure output.
76	!
77	! 667 2010-12-23 12:06:00Z suehring/gryschka
78	! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
79	! allocation of arrays. Calls of exchange_horiz are modified.
80	! Skip-value skip_do_avs changed to a dynamic adaption of ghost points.
81	!
82	! 646 2010-12-15 13:03:52Z raasch
83	! bugfix: missing define statements for netcdf added
84	!
85	! 493 2010-03-01 08:30:24Z raasch
86	! netCDF4 support (parallel output)
87	!
88	! 355 2009-07-17 01:03:01Z letzel
89	! simulated_time in netCDF output replaced by time_since_reference_point.
90	! Output of netCDF messages with aid of message handling routine.
91	! Output of messages replaced by message handling routine.
92	! Bugfix: to_be_resorted => s_av for time-averaged scalars
93	!
94	! 96 2007-06-04 08:07:41Z raasch
95	! Output of density and salinity
96	!
97	! 75 2007-03-22 09:54:05Z raasch
98	! 2nd+3rd argument removed from exchange horiz
99	!
100	! RCS Log replace by Id keyword, revision history cleaned up
101	!
102	! Revision 1.3 2006/06/02 15:18:59 raasch
103	! +argument "found", -argument grid in call of routine user_data_output_3d
104	!
105	! Revision 1.2 2006/02/23 10:23:07 raasch
106	! Former subroutine plot_3d renamed data_output_3d, pl.. renamed do..,
107	! .._anz renamed .._n,
108	! output extended to (almost) all quantities, output of user-defined quantities
109	!
110	! Revision 1.1 1997/09/03 06:29:36 raasch
111	! Initial revision
112	!
113	!
114	! Description:
115	! ------------
116	! Output of the 3D-arrays in netCDF and/or AVS format.
117	!------------------------------------------------------------------------------!
118
119	USE arrays_3d
120	USE averaging
121	USE cloud_parameters
122	USE control_parameters
123	USE cpulog
124	USE indices
125	USE interfaces
126	USE netcdf_control
127	USE particle_attributes
128	USE pegrid
129	USE precision_kind
130
131	IMPLICIT NONE
132
133	CHARACTER (LEN=9) :: simulated_time_mod
134
135	INTEGER :: av, i, if, j, k, n, pos, prec, psi
136
137	LOGICAL :: found, resorted
138
139	REAL :: mean_r, s_r3, s_r4
140
141	REAL(spk), DIMENSION(:,:,:), ALLOCATABLE :: local_pf
142
143	REAL, DIMENSION(:,:,:), POINTER :: to_be_resorted
144
145	!
146	!-- Return, if nothing to output
147	IF ( do3d_no(av) == 0 ) RETURN
148	!
149	!-- For parallel netcdf output the time axis must be limited. Return, if this
150	!-- limit is exceeded. This could be the case, if the simulated time exceeds
151	!-- the given end time by the length of the given output interval.
152	IF ( netcdf_data_format > 4 ) THEN
153	IF ( do3d_time_count(av) + 1 > ntdim_3d(av) ) THEN
154	WRITE ( message_string, * ) 'Output of 3d data is not given at t=', &
155	simulated_time, '&because the maximum ', &
156	'number of output time levels is ', &
157	'exceeded.'
158	CALL message( 'data_output_3d', 'PA0387', 0, 1, 0, 6, 0 )
159	RETURN
160	ENDIF
161	ENDIF
162
163	CALL cpu_log (log_point(14),'data_output_3d','start')
164
165	!
166	!-- Open output file.
167	!-- Also creates coordinate and fld-file for AVS.
168	!-- For classic or 64bit netCDF output or output of other (old) data formats,
169	!-- for a run on more than one PE, each PE opens its own file and
170	!-- writes the data of its subdomain in binary format (regardless of the format
171	!-- the user has requested). After the run, these files are combined to one
172	!-- file by combine_plot_fields in the format requested by the user (netcdf
173	!-- and/or avs).
174	!-- For netCDF4/HDF5 output, data is written in parallel into one file.
175	IF ( netcdf_output ) THEN
176	IF ( netcdf_data_format < 5 ) THEN
177	CALL check_open( 30 )
178	IF ( myid == 0 ) CALL check_open( 106+av*10 )
179	ELSE
180	CALL check_open( 106+av*10 )
181	ENDIF
182	ELSE
183	IF ( avs_output .OR. ( numprocs > 1 ) ) CALL check_open( 30 )
184	ENDIF
185
186	!
187	!-- Allocate a temporary array with the desired output dimensions.
188	ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb:nz_do3d) )
189
190	!
191	!-- Update the netCDF time axis
192	!-- In case of parallel output, this is only done by PE0 to increase the
193	!-- performance.
194	#if defined( __netcdf )
195	do3d_time_count(av) = do3d_time_count(av) + 1
196	IF ( myid == 0 ) THEN
197	IF ( netcdf_output ) THEN
198	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_time_3d(av), &
199	(/ time_since_reference_point /), &
200	start = (/ do3d_time_count(av) /), &
201	count = (/ 1 /) )
202	CALL handle_netcdf_error( 'data_output_3d', 376 )
203	ENDIF
204	ENDIF
205	#endif
206
207	!
208	!-- Loop over all variables to be written.
209	if = 1
210
211	DO WHILE ( do3d(av,if)(1:1) /= ' ' )
212	!
213	!-- Set the precision for data compression.
214	IF ( do3d_compress ) THEN
215	DO i = 1, 100
216	IF ( plot_3d_precision(i)%variable == do3d(av,if) ) THEN
217	prec = plot_3d_precision(i)%precision
218	EXIT
219	ENDIF
220	ENDDO
221	ENDIF
222
223	!
224	!-- Store the array chosen on the temporary array.
225	resorted = .FALSE.
226	SELECT CASE ( TRIM( do3d(av,if) ) )
227
228	CASE ( 'e' )
229	IF ( av == 0 ) THEN
230	to_be_resorted => e
231	ELSE
232	to_be_resorted => e_av
233	ENDIF
234
235	CASE ( 'lpt' )
236	IF ( av == 0 ) THEN
237	to_be_resorted => pt
238	ELSE
239	to_be_resorted => lpt_av
240	ENDIF
241
242	CASE ( 'nr' )
243	IF ( av == 0 ) THEN
244	to_be_resorted => nr
245	ELSE
246	to_be_resorted => nr_av
247	ENDIF
248
249	CASE ( 'p' )
250	IF ( av == 0 ) THEN
251	IF ( psolver /= 'sor' ) CALL exchange_horiz( p, nbgp )
252	to_be_resorted => p
253	ELSE
254	IF ( psolver /= 'sor' ) CALL exchange_horiz( p_av, nbgp )
255	to_be_resorted => p_av
256	ENDIF
257
258	CASE ( 'pc' ) ! particle concentration (requires ghostpoint exchange)
259	IF ( av == 0 ) THEN
260	tend = prt_count
261	CALL exchange_horiz( tend, nbgp )
262	DO i = nxlg, nxrg
263	DO j = nysg, nyng
264	DO k = nzb, nz_do3d
265	local_pf(i,j,k) = tend(k,j,i)
266	ENDDO
267	ENDDO
268	ENDDO
269	resorted = .TRUE.
270	ELSE
271	CALL exchange_horiz( pc_av, nbgp )
272	to_be_resorted => pc_av
273	ENDIF
274
275	CASE ( 'pr' ) ! mean particle radius
276	IF ( av == 0 ) THEN
277	DO i = nxl, nxr
278	DO j = nys, nyn
279	DO k = nzb, nz_do3d
280	psi = prt_start_index(k,j,i)
281	s_r3 = 0.0
282	s_r4 = 0.0
283	DO n = psi, psi+prt_count(k,j,i)-1
284	s_r3 = s_r3 + particles(n)%radius*3 &
285	particles(n)%weight_factor
286	s_r4 = s_r4 + particles(n)%radius*4 &
287	particles(n)%weight_factor
288	ENDDO
289	IF ( s_r3 /= 0.0 ) THEN
290	mean_r = s_r4 / s_r3
291	ELSE
292	mean_r = 0.0
293	ENDIF
294	tend(k,j,i) = mean_r
295	ENDDO
296	ENDDO
297	ENDDO
298	CALL exchange_horiz( tend, nbgp )
299	DO i = nxlg, nxrg
300	DO j = nysg, nyng
301	DO k = nzb, nz_do3d
302	local_pf(i,j,k) = tend(k,j,i)
303	ENDDO
304	ENDDO
305	ENDDO
306	resorted = .TRUE.
307	ELSE
308	CALL exchange_horiz( pr_av, nbgp )
309	to_be_resorted => pr_av
310	ENDIF
311
312	CASE ( 'prr' )
313	IF ( av == 0 ) THEN
314	CALL exchange_horiz( prr, nbgp )
315	DO i = nxlg, nxrg
316	DO j = nysg, nyng
317	DO k = nzb, nzt+1
318	local_pf(i,j,k) = prr(k,j,i)
319	ENDDO
320	ENDDO
321	ENDDO
322	ELSE
323	CALL exchange_horiz( prr_av, nbgp )
324	DO i = nxlg, nxrg
325	DO j = nysg, nyng
326	DO k = nzb, nzt+1
327	local_pf(i,j,k) = prr_av(k,j,i)
328	ENDDO
329	ENDDO
330	ENDDO
331	ENDIF
332	resorted = .TRUE.
333
334	CASE ( 'pt' )
335	IF ( av == 0 ) THEN
336	IF ( .NOT. cloud_physics ) THEN
337	to_be_resorted => pt
338	ELSE
339	DO i = nxlg, nxrg
340	DO j = nysg, nyng
341	DO k = nzb, nz_do3d
342	local_pf(i,j,k) = pt(k,j,i) + l_d_cp * &
343	pt_d_t(k) * &
344	ql(k,j,i)
345	ENDDO
346	ENDDO
347	ENDDO
348	resorted = .TRUE.
349	ENDIF
350	ELSE
351	to_be_resorted => pt_av
352	ENDIF
353
354	CASE ( 'q' )
355	IF ( av == 0 ) THEN
356	to_be_resorted => q
357	ELSE
358	to_be_resorted => q_av
359	ENDIF
360
361	CASE ( 'qc' )
362	IF ( av == 0 ) THEN
363	to_be_resorted => qc
364	ELSE
365	to_be_resorted => qc_av
366	ENDIF
367
368	CASE ( 'ql' )
369	IF ( av == 0 ) THEN
370	to_be_resorted => ql
371	ELSE
372	to_be_resorted => ql_av
373	ENDIF
374
375	CASE ( 'ql_c' )
376	IF ( av == 0 ) THEN
377	to_be_resorted => ql_c
378	ELSE
379	to_be_resorted => ql_c_av
380	ENDIF
381
382	CASE ( 'ql_v' )
383	IF ( av == 0 ) THEN
384	to_be_resorted => ql_v
385	ELSE
386	to_be_resorted => ql_v_av
387	ENDIF
388
389	CASE ( 'ql_vp' )
390	IF ( av == 0 ) THEN
391	DO i = nxl, nxr
392	DO j = nys, nyn
393	DO k = nzb, nz_do3d
394	psi = prt_start_index(k,j,i)
395	DO n = psi, psi+prt_count(k,j,i)-1
396	tend(k,j,i) = tend(k,j,i) + &
397	particles(n)%weight_factor / &
398	prt_count(k,j,i)
399	ENDDO
400	ENDDO
401	ENDDO
402	ENDDO
403	CALL exchange_horiz( tend, nbgp )
404	DO i = nxlg, nxrg
405	DO j = nysg, nyng
406	DO k = nzb, nz_do3d
407	local_pf(i,j,k) = tend(k,j,i)
408	ENDDO
409	ENDDO
410	ENDDO
411	resorted = .TRUE.
412	ELSE
413	CALL exchange_horiz( ql_vp_av, nbgp )
414	to_be_resorted => ql_vp_av
415	ENDIF
416
417	CASE ( 'qr' )
418	IF ( av == 0 ) THEN
419	to_be_resorted => qr
420	ELSE
421	to_be_resorted => qr_av
422	ENDIF
423
424	CASE ( 'qv' )
425	IF ( av == 0 ) THEN
426	DO i = nxlg, nxrg
427	DO j = nysg, nyng
428	DO k = nzb, nz_do3d
429	local_pf(i,j,k) = q(k,j,i) - ql(k,j,i)
430	ENDDO
431	ENDDO
432	ENDDO
433	resorted = .TRUE.
434	ELSE
435	to_be_resorted => qv_av
436	ENDIF
437
438	CASE ( 'rho' )
439	IF ( av == 0 ) THEN
440	to_be_resorted => rho
441	ELSE
442	to_be_resorted => rho_av
443	ENDIF
444
445	CASE ( 's' )
446	IF ( av == 0 ) THEN
447	to_be_resorted => q
448	ELSE
449	to_be_resorted => s_av
450	ENDIF
451
452	CASE ( 'sa' )
453	IF ( av == 0 ) THEN
454	to_be_resorted => sa
455	ELSE
456	to_be_resorted => sa_av
457	ENDIF
458
459	CASE ( 'u' )
460	IF ( av == 0 ) THEN
461	to_be_resorted => u
462	ELSE
463	to_be_resorted => u_av
464	ENDIF
465
466	CASE ( 'v' )
467	IF ( av == 0 ) THEN
468	to_be_resorted => v
469	ELSE
470	to_be_resorted => v_av
471	ENDIF
472
473	CASE ( 'vpt' )
474	IF ( av == 0 ) THEN
475	to_be_resorted => vpt
476	ELSE
477	to_be_resorted => vpt_av
478	ENDIF
479
480	CASE ( 'w' )
481	IF ( av == 0 ) THEN
482	to_be_resorted => w
483	ELSE
484	to_be_resorted => w_av
485	ENDIF
486
487	CASE DEFAULT
488	!
489	!-- User defined quantity
490	CALL user_data_output_3d( av, do3d(av,if), found, local_pf, &
491	nz_do3d )
492	resorted = .TRUE.
493
494	IF ( .NOT. found ) THEN
495	message_string = 'no output available for: ' // &
496	TRIM( do3d(av,if) )
497	CALL message( 'data_output_3d', 'PA0182', 0, 0, 0, 6, 0 )
498	ENDIF
499
500	END SELECT
501
502	!
503	!-- Resort the array to be output, if not done above
504	IF ( .NOT. resorted ) THEN
505	DO i = nxlg, nxrg
506	DO j = nysg, nyng
507	DO k = nzb, nz_do3d
508	local_pf(i,j,k) = to_be_resorted(k,j,i)
509	ENDDO
510	ENDDO
511	ENDDO
512	ENDIF
513
514	!
515	!-- Output of the volume data information for the AVS-FLD-file.
516	do3d_avs_n = do3d_avs_n + 1
517	IF ( myid == 0 .AND. avs_output ) THEN
518	!
519	!-- AVS-labels must not contain any colons. Hence they must be removed
520	!-- from the time character string.
521	simulated_time_mod = simulated_time_chr
522	DO WHILE ( SCAN( simulated_time_mod, ':' ) /= 0 )
523	pos = SCAN( simulated_time_mod, ':' )
524	simulated_time_mod(pos:pos) = '/'
525	ENDDO
526
527	IF ( av == 0 ) THEN
528	WRITE ( 33, 3300 ) do3d_avs_n, TRIM( avs_data_file ), &
529	skip_do_avs, TRIM( do3d(av,if) ), &
530	TRIM( simulated_time_mod )
531	ELSE
532	WRITE ( 33, 3300 ) do3d_avs_n, TRIM( avs_data_file ), &
533	skip_do_avs, TRIM( do3d(av,if) ) // &
534	' averaged', TRIM( simulated_time_mod )
535	ENDIF
536	!
537	!-- Determine the Skip-value for the next array. Record end and start
538	!-- require 4 byte each.
539	skip_do_avs = skip_do_avs + ( ( ( nx+2nbgp ) ( ny+2nbgp ) &
540	( nz_do3d+1 ) ) * 4 + 8 )
541	ENDIF
542
543	!
544	!-- Output of the 3D-array. (compressed/uncompressed)
545	IF ( do3d_compress ) THEN
546	!
547	!-- Compression, output of compression information on FLD-file and output
548	!-- of compressed data.
549	CALL write_compressed( local_pf, 30, 33, myid, nxl, nxr, nyn, nys, &
550	nzb, nz_do3d, prec, nbgp )
551	ELSE
552	!
553	!-- Uncompressed output.
554	#if defined( __parallel )
555	IF ( netcdf_output ) THEN
556	IF ( netcdf_data_format < 5 ) THEN
557	!
558	!-- Non-parallel netCDF output. Data is output in parallel in
559	!-- FORTRAN binary format here, and later collected into one file by
560	!-- combine_plot_fields
561	IF ( myid == 0 ) THEN
562	WRITE ( 30 ) time_since_reference_point, &
563	do3d_time_count(av), av
564	ENDIF
565	DO i = 0, io_blocks-1
566	IF ( i == io_group ) THEN
567	WRITE ( 30 ) nxlg, nxrg, nysg, nyng, nzb, nz_do3d
568	WRITE ( 30 ) local_pf
569	ENDIF
570	#if defined( __parallel )
571	CALL MPI_BARRIER( comm2d, ierr )
572	#endif
573	ENDDO
574
575	ELSE
576	#if defined( __netcdf )
577	!
578	!-- Parallel output in netCDF4/HDF5 format.
579	!-- Do not output redundant ghost point data except for the
580	!-- boundaries of the total domain.
581	IF ( nxr == nx .AND. nyn /= ny ) THEN
582	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if), &
583	local_pf(nxl:nxr+1,nys:nyn,nzb:nz_do3d), &
584	start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /), &
585	count = (/ nxr-nxl+2, nyn-nys+1, nz_do3d-nzb+1, 1 /) )
586	ELSEIF ( nxr /= nx .AND. nyn == ny ) THEN
587	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if), &
588	local_pf(nxl:nxr,nys:nyn+1,nzb:nz_do3d), &
589	start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /), &
590	count = (/ nxr-nxl+1, nyn-nys+2, nz_do3d-nzb+1, 1 /) )
591	ELSEIF ( nxr == nx .AND. nyn == ny ) THEN
592	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if), &
593	local_pf(nxl:nxr+1,nys:nyn+1,nzb:nz_do3d), &
594	start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /), &
595	count = (/ nxr-nxl+2, nyn-nys+2, nz_do3d-nzb+1, 1 /) )
596	ELSE
597	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if), &
598	local_pf(nxl:nxr,nys:nyn,nzb:nz_do3d), &
599	start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /), &
600	count = (/ nxr-nxl+1, nyn-nys+1, nz_do3d-nzb+1, 1 /) )
601	ENDIF
602	CALL handle_netcdf_error( 'data_output_3d', 386 )
603	#endif
604	ENDIF
605	ENDIF
606	#else
607	IF ( avs_output ) THEN
608	WRITE ( 30 ) local_pf(nxl:nxr+1,nys:nyn+1,:)
609	ENDIF
610	#if defined( __netcdf )
611	IF ( netcdf_output ) THEN
612
613	nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if), &
614	local_pf(nxl:nxr+1,nys:nyn+1,nzb:nz_do3d), &
615	start = (/ 1, 1, 1, do3d_time_count(av) /), &
616	count = (/ nx+2, ny+2, nz_do3d-nzb+1, 1 /) )
617	CALL handle_netcdf_error( 'data_output_3d', 446 )
618
619	ENDIF
620	#endif
621	#endif
622	ENDIF
623
624	if = if + 1
625
626	ENDDO
627
628	!
629	!-- Deallocate temporary array.
630	DEALLOCATE( local_pf )
631
632
633	CALL cpu_log (log_point(14),'data_output_3d','stop','nobarrier')
634
635	!
636	!-- Formats.
637	3300 FORMAT ('variable ',I4,' file=',A,' filetype=unformatted skip=',I12/ &
638	'label = ',A,A)
639
640	END SUBROUTINE data_output_3d

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |