Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |