Home

Context Navigation

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

Last change on this file since 1683 was 1683, checked in by knoop, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 22.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |