Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |