Home

Context Navigation

source: palm/trunk/SOURCE/data_output_mask.f90 @ 1438

Last change on this file since 1438 was 1438, checked in by heinze, 10 years ago
+nr, qc, qr for mask output
Property svn:keywords set to `Id`
File size: 21.5 KB

Line
1	SUBROUTINE data_output_mask( 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	! +nr, qc, qr
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: data_output_mask.f90 1438 2014-07-22 14:14:06Z heinze $
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	!
36	!
37	! 1320 2014-03-20 08:40:49Z raasch
38	! ONLY-attribute added to USE-statements,
39	! kind-parameters added to all INTEGER and REAL declaration statements,
40	! kinds are defined in new module kinds,
41	! revision history before 2012 removed,
42	! comment fields (!:) to be used for variable explanations added to
43	! all variable declaration statements
44	!
45	! 1318 2014-03-17 13:35:16Z raasch
46	! barrier argument removed from cpu_log,
47	! module interfaces removed
48	!
49	! 1092 2013-02-02 11:24:22Z raasch
50	! unused variables removed
51	!
52	! 1036 2012-10-22 13:43:42Z raasch
53	! code put under GPL (PALM 3.9)
54	!
55	! 1031 2012-10-19 14:35:30Z raasch
56	! netCDF4 without parallel file support implemented
57	!
58	! 1007 2012-09-19 14:30:36Z franke
59	! Bugfix: calculation of pr must depend on the particle weighting factor,
60	! missing calculation of ql_vp added
61	!
62	! 410 2009-12-04 17:05:40Z letzel
63	! Initial version
64	!
65	! Description:
66	! ------------
67	! Masked data output in netCDF format for current mask (current value of mid).
68	!------------------------------------------------------------------------------!
69
70	#if defined( __netcdf )
71	USE arrays_3d, &
72	ONLY: e, nr, p, pt, q, qc, ql, ql_c, ql_v, qr, rho, sa, tend, u, &
73	v, vpt, w
74
75	USE averaging, &
76	ONLY: e_av, lpt_av, nr_av, p_av, pc_av, pr_av, pt_av, q_av, qc_av, &
77	ql_av, ql_c_av, ql_v_av, ql_vp_av, qv_av, qr_av, rho_av, s_av, &
78	sa_av, u_av, v_av, vpt_av, w_av
79
80	USE cloud_parameters, &
81	ONLY: l_d_cp, pt_d_t
82
83	USE control_parameters, &
84	ONLY: cloud_physics, domask, domask_no, domask_time_count, &
85	mask_i, mask_j, mask_k, mask_size, mask_size_l, &
86	mask_start_l, max_masks, message_string, mid, &
87	netcdf_data_format, nz_do3d, simulated_time
88	USE cpulog, &
89	ONLY: cpu_log, log_point
90
91
92
93	USE indices, &
94	ONLY: nbgp, nxl, nxr, nyn, nys, nzb, nzt
95
96	USE kinds
97
98	USE netcdf
99
100	USE netcdf_control
101
102	USE particle_attributes, &
103	ONLY: grid_particles, number_of_particles, particles, &
104	particle_advection_start, prt_count
105
106	USE pegrid
107
108	IMPLICIT NONE
109
110	INTEGER(iwp) :: av !:
111	INTEGER(iwp) :: ngp !:
112	INTEGER(iwp) :: i !:
113	INTEGER(iwp) :: if !:
114	INTEGER(iwp) :: j !:
115	INTEGER(iwp) :: k !:
116	INTEGER(iwp) :: n !:
117	INTEGER(iwp) :: psi !:
118	INTEGER(iwp) :: sender !:
119	INTEGER(iwp) :: ind(6) !:
120
121	LOGICAL :: found !:
122	LOGICAL :: resorted !:
123
124	REAL(wp) :: mean_r !:
125	REAL(wp) :: s_r2 !:
126	REAL(wp) :: s_r3 !:
127
128	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: local_pf !:
129	#if defined( __parallel )
130	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: total_pf !:
131	#endif
132	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !:
133
134	!
135	!-- Return, if nothing to output
136	IF ( domask_no(mid,av) == 0 ) RETURN
137
138	CALL cpu_log (log_point(49),'data_output_mask','start')
139
140	!
141	!-- Open output file.
142	IF ( myid == 0 .OR. netcdf_data_format > 4 ) THEN
143	CALL check_open( 200+mid+av*max_masks )
144	ENDIF
145
146	!
147	!-- Allocate total and local output arrays.
148	#if defined( __parallel )
149	IF ( myid == 0 ) THEN
150	ALLOCATE( total_pf(mask_size(mid,1),mask_size(mid,2),mask_size(mid,3)) )
151	ENDIF
152	#endif
153	ALLOCATE( local_pf(mask_size_l(mid,1),mask_size_l(mid,2), &
154	mask_size_l(mid,3)) )
155
156	!
157	!-- Update the netCDF time axis.
158	domask_time_count(mid,av) = domask_time_count(mid,av) + 1
159	IF ( myid == 0 .OR. netcdf_data_format > 4 ) THEN
160	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_time_mask(mid,av), &
161	(/ simulated_time /), &
162	start = (/ domask_time_count(mid,av) /), &
163	count = (/ 1 /) )
164	CALL handle_netcdf_error( 'data_output_mask', 460 )
165	ENDIF
166
167	!
168	!-- Loop over all variables to be written.
169	if = 1
170
171	DO WHILE ( domask(mid,av,if)(1:1) /= ' ' )
172	!
173	!-- Reallocate local_pf on PE 0 since its shape changes during MPI exchange
174	IF ( netcdf_data_format < 5 .AND. myid == 0 .AND. if > 1 ) THEN
175	DEALLOCATE( local_pf )
176	ALLOCATE( local_pf(mask_size_l(mid,1),mask_size_l(mid,2), &
177	mask_size_l(mid,3)) )
178	ENDIF
179	!
180	!-- Store the variable chosen.
181	resorted = .FALSE.
182	SELECT CASE ( TRIM( domask(mid,av,if) ) )
183
184	CASE ( 'e' )
185	IF ( av == 0 ) THEN
186	to_be_resorted => e
187	ELSE
188	to_be_resorted => e_av
189	ENDIF
190
191	CASE ( 'lpt' )
192	IF ( av == 0 ) THEN
193	to_be_resorted => pt
194	ELSE
195	to_be_resorted => lpt_av
196	ENDIF
197
198	CASE ( 'nr' )
199	IF ( av == 0 ) THEN
200	to_be_resorted => nr
201	ELSE
202	to_be_resorted => nr_av
203	ENDIF
204
205	CASE ( 'p' )
206	IF ( av == 0 ) THEN
207	to_be_resorted => p
208	ELSE
209	to_be_resorted => p_av
210	ENDIF
211
212	CASE ( 'pc' ) ! particle concentration (requires ghostpoint exchange)
213	IF ( av == 0 ) THEN
214	tend = prt_count
215	CALL exchange_horiz( tend, nbgp )
216	DO i = 1, mask_size_l(mid,1)
217	DO j = 1, mask_size_l(mid,2)
218	DO k = 1, mask_size_l(mid,3)
219	local_pf(i,j,k) = tend(mask_k(mid,k), &
220	mask_j(mid,j),mask_i(mid,i))
221	ENDDO
222	ENDDO
223	ENDDO
224	resorted = .TRUE.
225	ELSE
226	CALL exchange_horiz( pc_av, nbgp )
227	to_be_resorted => pc_av
228	ENDIF
229
230	CASE ( 'pr' ) ! mean particle radius (effective radius)
231	IF ( av == 0 ) THEN
232	IF ( simulated_time >= particle_advection_start ) THEN
233	DO i = nxl, nxr
234	DO j = nys, nyn
235	DO k = nzb, nz_do3d
236	number_of_particles = prt_count(k,j,i)
237	IF (number_of_particles <= 0) CYCLE
238	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
239	s_r2 = 0.0_wp
240	s_r3 = 0.0_wp
241	DO n = 1, number_of_particles
242	IF ( particles(n)%particle_mask ) THEN
243	s_r2 = s_r2 + grid_particles(k,j,i)%particles(n)%radius*2 &
244	grid_particles(k,j,i)%particles(n)%weight_factor
245	s_r3 = s_r3 + grid_particles(k,j,i)%particles(n)%radius*3 &
246	grid_particles(k,j,i)%particles(n)%weight_factor
247	ENDIF
248	ENDDO
249	IF ( s_r2 > 0.0_wp ) THEN
250	mean_r = s_r3 / s_r2
251	ELSE
252	mean_r = 0.0_wp
253	ENDIF
254	tend(k,j,i) = mean_r
255	ENDDO
256	ENDDO
257	ENDDO
258	CALL exchange_horiz( tend, nbgp )
259	ELSE
260	tend = 0.0_wp
261	ENDIF
262	DO i = 1, mask_size_l(mid,1)
263	DO j = 1, mask_size_l(mid,2)
264	DO k = 1, mask_size_l(mid,3)
265	local_pf(i,j,k) = tend(mask_k(mid,k), &
266	mask_j(mid,j),mask_i(mid,i))
267	ENDDO
268	ENDDO
269	ENDDO
270	resorted = .TRUE.
271	ELSE
272	CALL exchange_horiz( pr_av, nbgp )
273	to_be_resorted => pr_av
274	ENDIF
275
276	CASE ( 'pt' )
277	IF ( av == 0 ) THEN
278	IF ( .NOT. cloud_physics ) THEN
279	to_be_resorted => pt
280	ELSE
281	DO i = 1, mask_size_l(mid,1)
282	DO j = 1, mask_size_l(mid,2)
283	DO k = 1, mask_size_l(mid,3)
284	local_pf(i,j,k) = &
285	pt(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i)) &
286	+ l_d_cp * pt_d_t(mask_k(mid,k)) * &
287	ql(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
288	ENDDO
289	ENDDO
290	ENDDO
291	resorted = .TRUE.
292	ENDIF
293	ELSE
294	to_be_resorted => pt_av
295	ENDIF
296
297	CASE ( 'q' )
298	IF ( av == 0 ) THEN
299	to_be_resorted => q
300	ELSE
301	to_be_resorted => q_av
302	ENDIF
303
304	CASE ( 'qc' )
305	IF ( av == 0 ) THEN
306	to_be_resorted => qc
307	ELSE
308	to_be_resorted => qc_av
309	ENDIF
310
311	CASE ( 'ql' )
312	IF ( av == 0 ) THEN
313	to_be_resorted => ql
314	ELSE
315	to_be_resorted => ql_av
316	ENDIF
317
318	CASE ( 'ql_c' )
319	IF ( av == 0 ) THEN
320	to_be_resorted => ql_c
321	ELSE
322	to_be_resorted => ql_c_av
323	ENDIF
324
325	CASE ( 'ql_v' )
326	IF ( av == 0 ) THEN
327	to_be_resorted => ql_v
328	ELSE
329	to_be_resorted => ql_v_av
330	ENDIF
331
332	CASE ( 'ql_vp' )
333	IF ( av == 0 ) THEN
334	IF ( simulated_time >= particle_advection_start ) THEN
335	DO i = nxl, nxr
336	DO j = nys, nyn
337	DO k = nzb, nz_do3d
338	number_of_particles = prt_count(k,j,i)
339	IF (number_of_particles <= 0) CYCLE
340	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
341	DO n = 1, number_of_particles
342	IF ( particles(n)%particle_mask ) THEN
343	tend(k,j,i) = tend(k,j,i) + &
344	particles(n)%weight_factor / &
345	prt_count(k,j,i)
346	ENDIF
347	ENDDO
348	ENDDO
349	ENDDO
350	ENDDO
351	CALL exchange_horiz( tend, nbgp )
352	ELSE
353	tend = 0.0_wp
354	ENDIF
355	DO i = 1, mask_size_l(mid,1)
356	DO j = 1, mask_size_l(mid,2)
357	DO k = 1, mask_size_l(mid,3)
358	local_pf(i,j,k) = tend(mask_k(mid,k), &
359	mask_j(mid,j),mask_i(mid,i))
360	ENDDO
361	ENDDO
362	ENDDO
363	resorted = .TRUE.
364	ELSE
365	CALL exchange_horiz( ql_vp_av, nbgp )
366	to_be_resorted => ql_vp_av
367	ENDIF
368
369	CASE ( 'qv' )
370	IF ( av == 0 ) THEN
371	DO i = 1, mask_size_l(mid,1)
372	DO j = 1, mask_size_l(mid,2)
373	DO k = 1, mask_size_l(mid,3)
374	local_pf(i,j,k) = &
375	q(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i)) - &
376	ql(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
377	ENDDO
378	ENDDO
379	ENDDO
380	resorted = .TRUE.
381	ELSE
382	to_be_resorted => qv_av
383	ENDIF
384
385	CASE ( 'qr' )
386	IF ( av == 0 ) THEN
387	to_be_resorted => qr
388	ELSE
389	to_be_resorted => qr_av
390	ENDIF
391
392	CASE ( 'rho' )
393	IF ( av == 0 ) THEN
394	to_be_resorted => rho
395	ELSE
396	to_be_resorted => rho_av
397	ENDIF
398
399	CASE ( 's' )
400	IF ( av == 0 ) THEN
401	to_be_resorted => q
402	ELSE
403	to_be_resorted => s_av
404	ENDIF
405
406	CASE ( 'sa' )
407	IF ( av == 0 ) THEN
408	to_be_resorted => sa
409	ELSE
410	to_be_resorted => sa_av
411	ENDIF
412
413	CASE ( 'u' )
414	IF ( av == 0 ) THEN
415	to_be_resorted => u
416	ELSE
417	to_be_resorted => u_av
418	ENDIF
419
420	CASE ( 'v' )
421	IF ( av == 0 ) THEN
422	to_be_resorted => v
423	ELSE
424	to_be_resorted => v_av
425	ENDIF
426
427	CASE ( 'vpt' )
428	IF ( av == 0 ) THEN
429	to_be_resorted => vpt
430	ELSE
431	to_be_resorted => vpt_av
432	ENDIF
433
434	CASE ( 'w' )
435	IF ( av == 0 ) THEN
436	to_be_resorted => w
437	ELSE
438	to_be_resorted => w_av
439	ENDIF
440
441	CASE DEFAULT
442	!
443	!-- User defined quantity
444	CALL user_data_output_mask(av, domask(mid,av,if), found, local_pf )
445	resorted = .TRUE.
446
447	IF ( .NOT. found ) THEN
448	WRITE ( message_string, * ) 'no output available for: ', &
449	TRIM( domask(mid,av,if) )
450	CALL message( 'data_output_mask', 'PA0327', 0, 0, 0, 6, 0 )
451	ENDIF
452
453	END SELECT
454
455	!
456	!-- Resort the array to be output, if not done above
457	IF ( .NOT. resorted ) THEN
458	DO i = 1, mask_size_l(mid,1)
459	DO j = 1, mask_size_l(mid,2)
460	DO k = 1, mask_size_l(mid,3)
461	local_pf(i,j,k) = to_be_resorted(mask_k(mid,k), &
462	mask_j(mid,j),mask_i(mid,i))
463	ENDDO
464	ENDDO
465	ENDDO
466	ENDIF
467
468	!
469	!-- I/O block. I/O methods are implemented
470	!-- (1) for parallel execution
471	!-- a. with netCDF 4 parallel I/O-enabled library
472	!-- b. with netCDF 3 library
473	!-- (2) for serial execution.
474	!-- The choice of method depends on the correct setting of preprocessor
475	!-- directives __parallel and __netcdf4_parallel as well as on the parameter
476	!-- netcdf_data_format.
477	#if defined( __parallel )
478	#if defined( __netcdf4_parallel )
479	IF ( netcdf_data_format > 4 ) THEN
480	!
481	!-- (1) a. Parallel I/O using netCDF 4 (not yet tested)
482	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
483	id_var_domask(mid,av,if), &
484	local_pf, &
485	start = (/ mask_start_l(mid,1), mask_start_l(mid,2), &
486	mask_start_l(mid,3), domask_time_count(mid,av) /), &
487	count = (/ mask_size_l(mid,1), mask_size_l(mid,2), &
488	mask_size_l(mid,3), 1 /) )
489	CALL handle_netcdf_error( 'data_output_mask', 461 )
490	ELSE
491	#endif
492	!
493	!-- (1) b. Conventional I/O only through PE0
494	!-- PE0 receives partial arrays from all processors of the respective mask
495	!-- and outputs them. Here a barrier has to be set, because otherwise
496	!-- "-MPI- FATAL: Remote protocol queue full" may occur.
497	CALL MPI_BARRIER( comm2d, ierr )
498
499	ngp = mask_size_l(mid,1) * mask_size_l(mid,2) * mask_size_l(mid,3)
500	IF ( myid == 0 ) THEN
501	!
502	!-- Local array can be relocated directly.
503	total_pf( &
504	mask_start_l(mid,1):mask_start_l(mid,1)+mask_size_l(mid,1)-1, &
505	mask_start_l(mid,2):mask_start_l(mid,2)+mask_size_l(mid,2)-1, &
506	mask_start_l(mid,3):mask_start_l(mid,3)+mask_size_l(mid,3)-1 ) &
507	= local_pf
508	!
509	!-- Receive data from all other PEs.
510	DO n = 1, numprocs-1
511	!
512	!-- Receive index limits first, then array.
513	!-- Index limits are received in arbitrary order from the PEs.
514	CALL MPI_RECV( ind(1), 6, MPI_INTEGER, MPI_ANY_SOURCE, 0, &
515	comm2d, status, ierr )
516	!
517	!-- Not all PEs have data for the mask
518	IF ( ind(1) /= -9999 ) THEN
519	ngp = ( ind(2)-ind(1)+1 ) * (ind(4)-ind(3)+1 ) * &
520	( ind(6)-ind(5)+1 )
521	sender = status(MPI_SOURCE)
522	DEALLOCATE( local_pf )
523	ALLOCATE(local_pf(ind(1):ind(2),ind(3):ind(4),ind(5):ind(6)))
524	CALL MPI_RECV( local_pf(ind(1),ind(3),ind(5)), ngp, &
525	MPI_REAL, sender, 1, comm2d, status, ierr )
526	total_pf(ind(1):ind(2),ind(3):ind(4),ind(5):ind(6)) &
527	= local_pf
528	ENDIF
529	ENDDO
530
531	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
532	id_var_domask(mid,av,if), total_pf, &
533	start = (/ 1, 1, 1, domask_time_count(mid,av) /), &
534	count = (/ mask_size(mid,1), mask_size(mid,2), &
535	mask_size(mid,3), 1 /) )
536	CALL handle_netcdf_error( 'data_output_mask', 462 )
537
538	ELSE
539	!
540	!-- If at least part of the mask resides on the PE, send the index
541	!-- limits for the target array, otherwise send -9999 to PE0.
542	IF ( mask_size_l(mid,1) > 0 .AND. mask_size_l(mid,2) > 0 .AND. &
543	mask_size_l(mid,3) > 0 ) &
544	THEN
545	ind(1) = mask_start_l(mid,1)
546	ind(2) = mask_start_l(mid,1) + mask_size_l(mid,1) - 1
547	ind(3) = mask_start_l(mid,2)
548	ind(4) = mask_start_l(mid,2) + mask_size_l(mid,2) - 1
549	ind(5) = mask_start_l(mid,3)
550	ind(6) = mask_start_l(mid,3) + mask_size_l(mid,3) - 1
551	ELSE
552	ind(1) = -9999; ind(2) = -9999
553	ind(3) = -9999; ind(4) = -9999
554	ind(5) = -9999; ind(6) = -9999
555	ENDIF
556	CALL MPI_SEND( ind(1), 6, MPI_INTEGER, 0, 0, comm2d, ierr )
557	!
558	!-- If applicable, send data to PE0.
559	IF ( ind(1) /= -9999 ) THEN
560	CALL MPI_SEND( local_pf(1,1,1), ngp, MPI_REAL, 0, 1, comm2d, &
561	ierr )
562	ENDIF
563	ENDIF
564	!
565	!-- A barrier has to be set, because otherwise some PEs may proceed too
566	!-- fast so that PE0 may receive wrong data on tag 0.
567	CALL MPI_BARRIER( comm2d, ierr )
568	#if defined( __netcdf4_parallel )
569	ENDIF
570	#endif
571	#else
572	!
573	!-- (2) For serial execution of PALM, the single processor (PE0) holds all
574	!-- data and writes them directly to file.
575	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
576	id_var_domask(mid,av,if), &
577	local_pf, &
578	start = (/ 1, 1, 1, domask_time_count(mid,av) /), &
579	count = (/ mask_size_l(mid,1), mask_size_l(mid,2), &
580	mask_size_l(mid,3), 1 /) )
581	CALL handle_netcdf_error( 'data_output_mask', 463 )
582	#endif
583
584	if = if + 1
585
586	ENDDO
587
588	!
589	!-- Deallocate temporary arrays.
590	DEALLOCATE( local_pf )
591	#if defined( __parallel )
592	IF ( myid == 0 ) THEN
593	DEALLOCATE( total_pf )
594	ENDIF
595	#endif
596
597
598	CALL cpu_log( log_point(49), 'data_output_mask', 'stop' )
599	#endif
600
601	END SUBROUTINE data_output_mask

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |