Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |