Home

Context Navigation

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

Last change on this file since 4256 was 4246, checked in by pavelkrc, 5 years ago
Convert CRLF line endings to LF
Property svn:keywords set to `Id`
File size: 30.8 KB

Line
1	!> @file data_output_mask.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! 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-2019 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: data_output_mask.f90 4246 2019-09-30 09:27:52Z monakurppa $
27	! Corrected "Former revisions" section
28	!
29	! 4168 2019-08-16 13:50:17Z suehring
30	! Remove variable grid
31	!
32	! 4167 2019-08-16 11:01:48Z suehring
33	! Changed behaviour of masked output over surface to follow terrain and ignore
34	! buildings (J.Resler, T.Gronemeier)
35	!
36	! 4069 2019-07-01 14:05:51Z Giersch
37	! Masked output running index mid has been introduced as a local variable to
38	! avoid runtime error (Loop variable has been modified) in time_integration
39	!
40	! 4039 2019-06-18 10:32:41Z suehring
41	! Modularize diagnostic output
42	!
43	! 3994 2019-05-22 18:08:09Z suehring
44	! output of turbulence intensity added
45	!
46	! 3665 2019-01-10 08:28:24Z raasch
47	! unused variables removed
48	!
49	! 3655 2019-01-07 16:51:22Z knoop
50	! Fix output time levels (use time_since_reference_point)
51	!
52	! 410 2009-12-04 17:05:40Z letzel
53	! Initial version
54	!
55	! Description:
56	! ------------
57	!> Masked data output in netCDF format for current mask (current value of mid).
58	!------------------------------------------------------------------------------!
59	SUBROUTINE data_output_mask( av, mid )
60
61
62
63	#if defined( __netcdf )
64	USE arrays_3d, &
65	ONLY: e, nc, nr, p, pt, q, qc, ql, ql_c, ql_v, qr, rho_ocean, s, sa, &
66	tend, u, v, vpt, w, d_exner
67
68	USE averaging, &
69	ONLY: e_av, lpt_av, nc_av, nr_av, p_av, pc_av, pr_av, pt_av, q_av, &
70	qc_av, ql_av, ql_c_av, ql_v_av, ql_vp_av, qv_av, qr_av, &
71	rho_ocean_av, s_av, sa_av, u_av, v_av, vpt_av, w_av
72
73	USE basic_constants_and_equations_mod, &
74	ONLY: lv_d_cp
75
76	USE chemistry_model_mod, &
77	ONLY: chem_data_output_mask
78
79	USE control_parameters, &
80	ONLY: air_chemistry, domask, domask_no, domask_time_count, mask_i, &
81	mask_j, mask_k, mask_size, mask_size_l, mask_start_l, &
82	mask_surface, &
83	max_masks, message_string, nz_do3d, salsa, &
84	time_since_reference_point
85
86	USE diagnostic_output_quantities_mod, &
87	ONLY: doq_output_mask
88
89	USE cpulog, &
90	ONLY: cpu_log, log_point
91
92	USE indices, &
93	ONLY: nbgp, nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
94
95	USE kinds
96
97	USE bulk_cloud_model_mod, &
98	ONLY: bulk_cloud_model
99
100	USE NETCDF
101
102	USE netcdf_interface, &
103	ONLY: fill_value, id_set_mask, id_var_domask, id_var_time_mask, &
104	nc_stat, netcdf_data_format, netcdf_handle_error
105
106	USE particle_attributes, &
107	ONLY: grid_particles, number_of_particles, particles, &
108	particle_advection_start, prt_count
109
110	USE pegrid
111
112	USE radiation_model_mod, &
113	ONLY: radiation, radiation_data_output_mask
114
115	USE salsa_mod, &
116	ONLY: salsa_data_output_mask
117
118
119	IMPLICIT NONE
120
121	INTEGER(iwp) :: av !< flag for (non-)average output
122	INTEGER(iwp) :: ngp !< number of grid points of an output slice
123	INTEGER(iwp) :: i !< loop index
124	INTEGER(iwp) :: ivar !< variable index
125	INTEGER(iwp) :: j !< loop index
126	INTEGER(iwp) :: k !< loop index
127	INTEGER(iwp) :: im !< loop index for masked variables
128	INTEGER(iwp) :: jm !< loop index for masked variables
129	INTEGER(iwp) :: kk !< vertical index
130	INTEGER(iwp) :: mid !< masked output running index
131	INTEGER(iwp) :: n !< loop index
132	INTEGER(iwp) :: netcdf_data_format_save !< value of netcdf_data_format
133	INTEGER(iwp) :: sender !< PE id of sending PE
134	INTEGER(iwp) :: ktt !< k index of highest terrain surface
135	INTEGER(iwp) :: ind(6) !< index limits (lower/upper bounds) of array 'local_2d'
136
137	LOGICAL :: found !< true if output variable was found
138	LOGICAL :: resorted !< true if variable is resorted
139
140	REAL(wp) :: mean_r !< mean particle radius
141	REAL(wp) :: s_r2 !< sum( particle-radius**2 )
142	REAL(wp) :: s_r3 !< sum( particle-radius**3 )
143
144	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: local_pf !< output array
145	#if defined( __parallel )
146	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: total_pf !< collected output array
147	#endif
148	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which shall be output
149
150	!
151	!-- Return, if nothing to output
152	IF ( domask_no(mid,av) == 0 ) RETURN
153
154	CALL cpu_log (log_point(49),'data_output_mask','start')
155
156	!
157	!-- Parallel netcdf output is not tested so far for masked data, hence
158	!-- netcdf_data_format is switched back to non-paralell output.
159	netcdf_data_format_save = netcdf_data_format
160	IF ( netcdf_data_format == 5 ) netcdf_data_format = 3
161	IF ( netcdf_data_format == 6 ) netcdf_data_format = 4
162
163	!
164	!-- Open output file.
165	IF ( myid == 0 .OR. netcdf_data_format > 4 ) THEN
166	CALL check_open( 200+mid+av*max_masks )
167	ENDIF
168
169	!
170	!-- Allocate total and local output arrays.
171	#if defined( __parallel )
172	IF ( myid == 0 ) THEN
173	ALLOCATE( total_pf(mask_size(mid,1),mask_size(mid,2),mask_size(mid,3)) )
174	ENDIF
175	#endif
176	ALLOCATE( local_pf(mask_size_l(mid,1),mask_size_l(mid,2), &
177	mask_size_l(mid,3)) )
178
179	!
180	!-- Update the netCDF time axis.
181	domask_time_count(mid,av) = domask_time_count(mid,av) + 1
182	IF ( myid == 0 .OR. netcdf_data_format > 4 ) THEN
183	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_time_mask(mid,av), &
184	(/ time_since_reference_point /), &
185	start = (/ domask_time_count(mid,av) /), &
186	count = (/ 1 /) )
187	CALL netcdf_handle_error( 'data_output_mask', 460 )
188	ENDIF
189
190	!
191	!-- Loop over all variables to be written.
192	ivar = 1
193
194	DO WHILE ( domask(mid,av,ivar)(1:1) /= ' ' )
195	!
196	!-- Reallocate local_pf on PE 0 since its shape changes during MPI exchange
197	IF ( netcdf_data_format < 5 .AND. myid == 0 .AND. ivar > 1 ) THEN
198	DEALLOCATE( local_pf )
199	ALLOCATE( local_pf(mask_size_l(mid,1),mask_size_l(mid,2), &
200	mask_size_l(mid,3)) )
201	ENDIF
202	!
203	!-- Store the variable chosen.
204	resorted = .FALSE.
205	SELECT CASE ( TRIM( domask(mid,av,ivar) ) )
206
207	CASE ( 'e' )
208	IF ( av == 0 ) THEN
209	to_be_resorted => e
210	ELSE
211	to_be_resorted => e_av
212	ENDIF
213
214	CASE ( 'thetal' )
215	IF ( av == 0 ) THEN
216	to_be_resorted => pt
217	ELSE
218	to_be_resorted => lpt_av
219	ENDIF
220
221	CASE ( 'nc' )
222	IF ( av == 0 ) THEN
223	to_be_resorted => nc
224	ELSE
225	to_be_resorted => nc_av
226	ENDIF
227
228	CASE ( 'nr' )
229	IF ( av == 0 ) THEN
230	to_be_resorted => nr
231	ELSE
232	to_be_resorted => nr_av
233	ENDIF
234
235	CASE ( 'p' )
236	IF ( av == 0 ) THEN
237	to_be_resorted => p
238	ELSE
239	to_be_resorted => p_av
240	ENDIF
241
242	CASE ( 'pc' ) ! particle concentration (requires ghostpoint exchange)
243	IF ( av == 0 ) THEN
244	tend = prt_count
245	CALL exchange_horiz( tend, nbgp )
246	IF ( .NOT. mask_surface(mid) ) THEN
247	DO i = 1, mask_size_l(mid,1)
248	DO j = 1, mask_size_l(mid,2)
249	DO k = 1, mask_size_l(mid,3)
250	local_pf(i,j,k) = tend(mask_k(mid,k), &
251	mask_j(mid,j),mask_i(mid,i))
252	ENDDO
253	ENDDO
254	ENDDO
255	ELSE
256	!
257	!-- Terrain-following masked output
258	DO i = 1, mask_size_l(mid,1)
259	DO j = 1, mask_size_l(mid,2)
260	!-- Get k index of the highest terraing surface
261	im = mask_i(mid,i)
262	jm = mask_j(mid,j)
263	ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )),&
264	DIM = 1 ) - 1
265	DO k = 1, mask_size_l(mid,3)
266	kk = MIN( ktt+mask_k(mid,k), nzt+1 )
267	!-- Set value if not in building
268	IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) ) THEN
269	local_pf(i,j,k) = fill_value
270	ELSE
271	local_pf(i,j,k) = tend(kk,jm,im)
272	ENDIF
273	ENDDO
274	ENDDO
275	ENDDO
276	ENDIF
277	resorted = .TRUE.
278	ELSE
279	CALL exchange_horiz( pc_av, nbgp )
280	to_be_resorted => pc_av
281	ENDIF
282
283	CASE ( 'pr' ) ! mean particle radius (effective radius)
284	IF ( av == 0 ) THEN
285	IF ( time_since_reference_point >= particle_advection_start ) THEN
286	DO i = nxl, nxr
287	DO j = nys, nyn
288	DO k = nzb, nz_do3d
289	number_of_particles = prt_count(k,j,i)
290	IF (number_of_particles <= 0) CYCLE
291	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
292	s_r2 = 0.0_wp
293	s_r3 = 0.0_wp
294	DO n = 1, number_of_particles
295	IF ( particles(n)%particle_mask ) THEN
296	s_r2 = s_r2 + grid_particles(k,j,i)%particles(n)%radius*2 &
297	grid_particles(k,j,i)%particles(n)%weight_factor
298	s_r3 = s_r3 + grid_particles(k,j,i)%particles(n)%radius*3 &
299	grid_particles(k,j,i)%particles(n)%weight_factor
300	ENDIF
301	ENDDO
302	IF ( s_r2 > 0.0_wp ) THEN
303	mean_r = s_r3 / s_r2
304	ELSE
305	mean_r = 0.0_wp
306	ENDIF
307	tend(k,j,i) = mean_r
308	ENDDO
309	ENDDO
310	ENDDO
311	CALL exchange_horiz( tend, nbgp )
312	ELSE
313	tend = 0.0_wp
314	ENDIF
315	IF ( .NOT. mask_surface(mid) ) THEN
316	DO i = 1, mask_size_l(mid,1)
317	DO j = 1, mask_size_l(mid,2)
318	DO k = 1, mask_size_l(mid,3)
319	local_pf(i,j,k) = tend(mask_k(mid,k), &
320	mask_j(mid,j),mask_i(mid,i))
321	ENDDO
322	ENDDO
323	ENDDO
324	ELSE
325	!
326	!-- Terrain-following masked output
327	DO i = 1, mask_size_l(mid,1)
328	DO j = 1, mask_size_l(mid,2)
329	!-- Get k index of the highest terraing surface
330	im = mask_i(mid,i)
331	jm = mask_j(mid,j)
332	ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
333	DO k = 1, mask_size_l(mid,3)
334	kk = MIN( ktt+mask_k(mid,k), nzt+1 )
335	!-- Set value if not in building
336	IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) ) THEN
337	local_pf(i,j,k) = fill_value
338	ELSE
339	local_pf(i,j,k) = tend(kk,jm,im)
340	ENDIF
341	ENDDO
342	ENDDO
343	ENDDO
344	ENDIF
345	resorted = .TRUE.
346	ELSE
347	CALL exchange_horiz( pr_av, nbgp )
348	to_be_resorted => pr_av
349	ENDIF
350
351	CASE ( 'theta' )
352	IF ( av == 0 ) THEN
353	IF ( .NOT. bulk_cloud_model ) THEN
354	to_be_resorted => pt
355	ELSE
356	IF ( .NOT. mask_surface(mid) ) THEN
357	DO i = 1, mask_size_l(mid,1)
358	DO j = 1, mask_size_l(mid,2)
359	DO k = 1, mask_size_l(mid,3)
360	local_pf(i,j,k) = &
361	pt(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i)) &
362	+ lv_d_cp * d_exner(mask_k(mid,k)) * &
363	ql(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
364	ENDDO
365	ENDDO
366	ENDDO
367	ELSE
368	!
369	!-- Terrain-following masked output
370	DO i = 1, mask_size_l(mid,1)
371	DO j = 1, mask_size_l(mid,2)
372	!-- Get k index of the highest terraing surface
373	im = mask_i(mid,i)
374	jm = mask_j(mid,j)
375	ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
376	DO k = 1, mask_size_l(mid,3)
377	kk = MIN( ktt+mask_k(mid,k), nzt+1 )
378	!-- Set value if not in building
379	IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) ) THEN
380	local_pf(i,j,k) = fill_value
381	ELSE
382	local_pf(i,j,k) = pt(kk,jm,im) + lv_d_cp * d_exner(kk) * ql(kk,jm,im)
383	ENDIF
384	ENDDO
385	ENDDO
386	ENDDO
387	ENDIF
388	resorted = .TRUE.
389	ENDIF
390	ELSE
391	to_be_resorted => pt_av
392	ENDIF
393
394	CASE ( 'q' )
395	IF ( av == 0 ) THEN
396	to_be_resorted => q
397	ELSE
398	to_be_resorted => q_av
399	ENDIF
400
401	CASE ( 'qc' )
402	IF ( av == 0 ) THEN
403	to_be_resorted => qc
404	ELSE
405	to_be_resorted => qc_av
406	ENDIF
407
408	CASE ( 'ql' )
409	IF ( av == 0 ) THEN
410	to_be_resorted => ql
411	ELSE
412	to_be_resorted => ql_av
413	ENDIF
414
415	CASE ( 'ql_c' )
416	IF ( av == 0 ) THEN
417	to_be_resorted => ql_c
418	ELSE
419	to_be_resorted => ql_c_av
420	ENDIF
421
422	CASE ( 'ql_v' )
423	IF ( av == 0 ) THEN
424	to_be_resorted => ql_v
425	ELSE
426	to_be_resorted => ql_v_av
427	ENDIF
428
429	CASE ( 'ql_vp' )
430	IF ( av == 0 ) THEN
431	IF ( time_since_reference_point >= particle_advection_start ) THEN
432	DO i = nxl, nxr
433	DO j = nys, nyn
434	DO k = nzb, nz_do3d
435	number_of_particles = prt_count(k,j,i)
436	IF (number_of_particles <= 0) CYCLE
437	particles => grid_particles(k,j,i)%particles(1:number_of_particles)
438	DO n = 1, number_of_particles
439	IF ( particles(n)%particle_mask ) THEN
440	tend(k,j,i) = tend(k,j,i) + &
441	particles(n)%weight_factor / &
442	prt_count(k,j,i)
443	ENDIF
444	ENDDO
445	ENDDO
446	ENDDO
447	ENDDO
448	CALL exchange_horiz( tend, nbgp )
449	ELSE
450	tend = 0.0_wp
451	ENDIF
452	IF ( .NOT. mask_surface(mid) ) THEN
453	DO i = 1, mask_size_l(mid,1)
454	DO j = 1, mask_size_l(mid,2)
455	DO k = 1, mask_size_l(mid,3)
456	local_pf(i,j,k) = tend(mask_k(mid,k), &
457	mask_j(mid,j),mask_i(mid,i))
458	ENDDO
459	ENDDO
460	ENDDO
461	ELSE
462	!
463	!-- Terrain-following masked output
464	DO i = 1, mask_size_l(mid,1)
465	DO j = 1, mask_size_l(mid,2)
466	!-- Get k index of the highest terraing surface
467	im = mask_i(mid,i)
468	jm = mask_j(mid,j)
469	ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
470	DO k = 1, mask_size_l(mid,3)
471	kk = MIN( ktt+mask_k(mid,k), nzt+1 )
472	!-- Set value if not in building
473	IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) ) THEN
474	local_pf(i,j,k) = fill_value
475	ELSE
476	local_pf(i,j,k) = tend(kk,jm,im)
477	ENDIF
478	ENDDO
479	ENDDO
480	ENDDO
481	ENDIF
482	resorted = .TRUE.
483	ELSE
484	CALL exchange_horiz( ql_vp_av, nbgp )
485	to_be_resorted => ql_vp_av
486	ENDIF
487
488	CASE ( 'qv' )
489	IF ( av == 0 ) THEN
490	IF ( .NOT. mask_surface(mid) ) THEN
491	DO i = 1, mask_size_l(mid,1)
492	DO j = 1, mask_size_l(mid,2)
493	DO k = 1, mask_size_l(mid,3)
494	local_pf(i,j,k) = &
495	q(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i)) - &
496	ql(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
497	ENDDO
498	ENDDO
499	ENDDO
500	ELSE
501	!
502	!-- Terrain-following masked output
503	DO i = 1, mask_size_l(mid,1)
504	DO j = 1, mask_size_l(mid,2)
505	!-- Get k index of the highest terraing surface
506	im = mask_i(mid,i)
507	jm = mask_j(mid,j)
508	ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
509	DO k = 1, mask_size_l(mid,3)
510	kk = MIN( ktt+mask_k(mid,k), nzt+1 )
511	!-- Set value if not in building
512	IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) ) THEN
513	local_pf(i,j,k) = fill_value
514	ELSE
515	local_pf(i,j,k) = q(kk,jm,im) - ql(kk,jm,im)
516	ENDIF
517	ENDDO
518	ENDDO
519	ENDDO
520	ENDIF
521	resorted = .TRUE.
522	ELSE
523	to_be_resorted => qv_av
524	ENDIF
525
526	CASE ( 'qr' )
527	IF ( av == 0 ) THEN
528	to_be_resorted => qr
529	ELSE
530	to_be_resorted => qr_av
531	ENDIF
532
533	CASE ( 'rho_sea_water' )
534	IF ( av == 0 ) THEN
535	to_be_resorted => rho_ocean
536	ELSE
537	to_be_resorted => rho_ocean_av
538	ENDIF
539
540	CASE ( 's' )
541	IF ( av == 0 ) THEN
542	to_be_resorted => s
543	ELSE
544	to_be_resorted => s_av
545	ENDIF
546
547	CASE ( 'sa' )
548	IF ( av == 0 ) THEN
549	to_be_resorted => sa
550	ELSE
551	to_be_resorted => sa_av
552	ENDIF
553
554	CASE ( 'u' )
555	IF ( av == 0 ) THEN
556	to_be_resorted => u
557	ELSE
558	to_be_resorted => u_av
559	ENDIF
560
561	CASE ( 'v' )
562	IF ( av == 0 ) THEN
563	to_be_resorted => v
564	ELSE
565	to_be_resorted => v_av
566	ENDIF
567
568	CASE ( 'thetav' )
569	IF ( av == 0 ) THEN
570	to_be_resorted => vpt
571	ELSE
572	to_be_resorted => vpt_av
573	ENDIF
574
575	CASE ( 'w' )
576	IF ( av == 0 ) THEN
577	to_be_resorted => w
578	ELSE
579	to_be_resorted => w_av
580	ENDIF
581
582	CASE DEFAULT
583	!
584	!-- Set flag to steer output of radiation, land-surface, or user-defined
585	!-- quantities
586	found = .FALSE.
587	!
588	!-- Radiation quantity
589	IF ( .NOT. found .AND. radiation ) THEN
590	CALL radiation_data_output_mask(av, domask(mid,av,ivar), found,&
591	local_pf, mid )
592	ENDIF
593
594	IF ( .NOT. found .AND. air_chemistry ) THEN
595	CALL chem_data_output_mask(av, domask(mid,av,ivar), found, &
596	local_pf, mid )
597	ENDIF
598	!
599	!-- Check for diagnostic quantities
600	IF ( .NOT. found ) THEN
601	CALL doq_output_mask( av, domask(mid,av,ivar), found, local_pf, &
602	mid)
603	ENDIF
604	!
605	!-- SALSA quantities
606	IF ( .NOT. found .AND. salsa ) THEN
607	CALL salsa_data_output_mask( av, domask(mid,av,ivar), found, &
608	local_pf, mid )
609	ENDIF
610	!
611	!-- User defined quantity
612	IF ( .NOT. found ) THEN
613	CALL user_data_output_mask(av, domask(mid,av,ivar), found, &
614	local_pf, mid )
615	ENDIF
616
617	resorted = .TRUE.
618
619	IF ( .NOT. found ) THEN
620	WRITE ( message_string, * ) 'no masked output available for: ',&
621	TRIM( domask(mid,av,ivar) )
622	CALL message( 'data_output_mask', 'PA0327', 0, 0, 0, 6, 0 )
623	ENDIF
624
625	END SELECT
626
627	!
628	!-- Resort the array to be output, if not done above
629	IF ( .NOT. resorted ) THEN
630	IF ( .NOT. mask_surface(mid) ) THEN
631	!
632	!-- Default masked output
633	DO i = 1, mask_size_l(mid,1)
634	DO j = 1, mask_size_l(mid,2)
635	DO k = 1, mask_size_l(mid,3)
636	local_pf(i,j,k) = to_be_resorted(mask_k(mid,k), &
637	mask_j(mid,j),mask_i(mid,i))
638	ENDDO
639	ENDDO
640	ENDDO
641
642	ELSE
643	!
644	!-- Terrain-following masked output
645	DO i = 1, mask_size_l(mid,1)
646	DO j = 1, mask_size_l(mid,2)
647	!-- Get k index of the highest terraing surface
648	im = mask_i(mid,i)
649	jm = mask_j(mid,j)
650	ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_0(:,jm,im), 5 )), DIM = 1 ) - 1
651	DO k = 1, mask_size_l(mid,3)
652	kk = MIN( ktt+mask_k(mid,k), nzt+1 )
653	!-- Set value if not in building
654	IF ( BTEST( wall_flags_0(kk,jm,im), 6 ) ) THEN
655	local_pf(i,j,k) = fill_value
656	ELSE
657	local_pf(i,j,k) = to_be_resorted(kk,jm,im)
658	ENDIF
659	ENDDO
660	ENDDO
661	ENDDO
662
663	ENDIF
664	ENDIF
665
666	!
667	!-- I/O block. I/O methods are implemented
668	!-- (1) for parallel execution
669	!-- a. with netCDF 4 parallel I/O-enabled library
670	!-- b. with netCDF 3 library
671	!-- (2) for serial execution.
672	!-- The choice of method depends on the correct setting of preprocessor
673	!-- directives __parallel and __netcdf4_parallel as well as on the parameter
674	!-- netcdf_data_format.
675	#if defined( __parallel )
676	#if defined( __netcdf4_parallel )
677	IF ( netcdf_data_format > 4 ) THEN
678	!
679	!-- (1) a. Parallel I/O using netCDF 4 (not yet tested)
680	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
681	id_var_domask(mid,av,ivar), local_pf, &
682	start = (/ mask_start_l(mid,1), mask_start_l(mid,2), &
683	mask_start_l(mid,3), domask_time_count(mid,av) /), &
684	count = (/ mask_size_l(mid,1), mask_size_l(mid,2), &
685	mask_size_l(mid,3), 1 /) )
686	CALL netcdf_handle_error( 'data_output_mask', 461 )
687	ELSE
688	#endif
689	!
690	!-- (1) b. Conventional I/O only through PE0
691	!-- PE0 receives partial arrays from all processors of the respective mask
692	!-- and outputs them. Here a barrier has to be set, because otherwise
693	!-- "-MPI- FATAL: Remote protocol queue full" may occur.
694	CALL MPI_BARRIER( comm2d, ierr )
695
696	ngp = mask_size_l(mid,1) * mask_size_l(mid,2) * mask_size_l(mid,3)
697	IF ( myid == 0 ) THEN
698	!
699	!-- Local array can be relocated directly.
700	total_pf( &
701	mask_start_l(mid,1):mask_start_l(mid,1)+mask_size_l(mid,1)-1, &
702	mask_start_l(mid,2):mask_start_l(mid,2)+mask_size_l(mid,2)-1, &
703	mask_start_l(mid,3):mask_start_l(mid,3)+mask_size_l(mid,3)-1 ) &
704	= local_pf
705	!
706	!-- Receive data from all other PEs.
707	DO n = 1, numprocs-1
708	!
709	!-- Receive index limits first, then array.
710	!-- Index limits are received in arbitrary order from the PEs.
711	CALL MPI_RECV( ind(1), 6, MPI_INTEGER, MPI_ANY_SOURCE, 0, &
712	comm2d, status, ierr )
713	!
714	!-- Not all PEs have data for the mask
715	IF ( ind(1) /= -9999 ) THEN
716	ngp = ( ind(2)-ind(1)+1 ) * (ind(4)-ind(3)+1 ) * &
717	( ind(6)-ind(5)+1 )
718	sender = status(MPI_SOURCE)
719	DEALLOCATE( local_pf )
720	ALLOCATE(local_pf(ind(1):ind(2),ind(3):ind(4),ind(5):ind(6)))
721	CALL MPI_RECV( local_pf(ind(1),ind(3),ind(5)), ngp, &
722	MPI_REAL, sender, 1, comm2d, status, ierr )
723	total_pf(ind(1):ind(2),ind(3):ind(4),ind(5):ind(6)) &
724	= local_pf
725	ENDIF
726	ENDDO
727
728	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
729	id_var_domask(mid,av,ivar), total_pf, &
730	start = (/ 1, 1, 1, domask_time_count(mid,av) /), &
731	count = (/ mask_size(mid,1), mask_size(mid,2), &
732	mask_size(mid,3), 1 /) )
733	CALL netcdf_handle_error( 'data_output_mask', 462 )
734
735	ELSE
736	!
737	!-- If at least part of the mask resides on the PE, send the index
738	!-- limits for the target array, otherwise send -9999 to PE0.
739	IF ( mask_size_l(mid,1) > 0 .AND. mask_size_l(mid,2) > 0 .AND. &
740	mask_size_l(mid,3) > 0 ) &
741	THEN
742	ind(1) = mask_start_l(mid,1)
743	ind(2) = mask_start_l(mid,1) + mask_size_l(mid,1) - 1
744	ind(3) = mask_start_l(mid,2)
745	ind(4) = mask_start_l(mid,2) + mask_size_l(mid,2) - 1
746	ind(5) = mask_start_l(mid,3)
747	ind(6) = mask_start_l(mid,3) + mask_size_l(mid,3) - 1
748	ELSE
749	ind(1) = -9999; ind(2) = -9999
750	ind(3) = -9999; ind(4) = -9999
751	ind(5) = -9999; ind(6) = -9999
752	ENDIF
753	CALL MPI_SEND( ind(1), 6, MPI_INTEGER, 0, 0, comm2d, ierr )
754	!
755	!-- If applicable, send data to PE0.
756	IF ( ind(1) /= -9999 ) THEN
757	CALL MPI_SEND( local_pf(1,1,1), ngp, MPI_REAL, 0, 1, comm2d, &
758	ierr )
759	ENDIF
760	ENDIF
761	!
762	!-- A barrier has to be set, because otherwise some PEs may proceed too
763	!-- fast so that PE0 may receive wrong data on tag 0.
764	CALL MPI_BARRIER( comm2d, ierr )
765	#if defined( __netcdf4_parallel )
766	ENDIF
767	#endif
768	#else
769	!
770	!-- (2) For serial execution of PALM, the single processor (PE0) holds all
771	!-- data and writes them directly to file.
772	nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), &
773	id_var_domask(mid,av,ivar), local_pf, &
774	start = (/ 1, 1, 1, domask_time_count(mid,av) /), &
775	count = (/ mask_size_l(mid,1), mask_size_l(mid,2),&
776	mask_size_l(mid,3), 1 /) )
777	CALL netcdf_handle_error( 'data_output_mask', 463 )
778	#endif
779
780	ivar = ivar + 1
781
782	ENDDO
783
784	!
785	!-- Deallocate temporary arrays.
786	DEALLOCATE( local_pf )
787	#if defined( __parallel )
788	IF ( myid == 0 ) THEN
789	DEALLOCATE( total_pf )
790	ENDIF
791	#endif
792
793	!
794	!-- Switch back to original format given by user (see beginning of this routine)
795	netcdf_data_format = netcdf_data_format_save
796
797	CALL cpu_log( log_point(49), 'data_output_mask', 'stop' )
798	#endif
799
800
801	END SUBROUTINE data_output_mask

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |