Home

Context Navigation

data_output_3d.f90 @ 1327

Last change on this file since 1327 was 1327, checked in by raasch, 11 years ago

Changed:

-s real64 removed (.mrun.config.hlrnIII)
-r8 removed (.mrun.config.imuk)
deleted: .mrun.config.imuk_ice2_netcdf4 .mrun.config.imuk_hlrn

REAL constants defined as wp-kind in modules

"baroclinicity" renamed "baroclinity", "ocean version" replaced by
"ocean mode"

code parts concerning old output formats "iso2d" and "avs" removed.
netCDF is the only remaining output format.

Errors:

bugfix: duplicate error message 56 removed

Property svn:keywords set to Id

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

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

Context Navigation

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

Download in other formats: