Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |