Home

Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

diagnostic_output_quantities_mod.f90

Last change on this file 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: 97.3 KB

Line
1	!> @file diagnostic_output_quantities_mod.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	! Copyright 1997-2021 Leibniz Universitaet Hannover
17	!--------------------------------------------------------------------------------------------------!
18	!
19	! Current revisions:
20	! ------------------
21	!
22	!
23	! Former revisions:
24	! -----------------
25	! $Id: diagnostic_output_quantities_mod.f90 4895 2021-03-03 15:39:08Z banzhafs $
26	! Remove offset in terrain-following masked output
27	!
28	! 4866 2021-02-05 17:50:14Z suehring
29	! Implemented vertical passive scalar flux
30	! Implemented air temperature as diagnostic output quantity (previous commit)
31	!
32	! 4861 2021-02-01 10:51:45Z raasch
33	! Enable 3D data output also with 64-bit precision
34	!
35	! 4757 2020-10-26 10:23:38Z schwenkel
36	! Implement relative humidity as diagnostic output quantity
37	!
38	! 4671 2020-09-09 20:27:58Z pavelkrc
39	! Implementation of downward facing USM and LSM surfaces
40	!
41	! 4583 2020-06-29 12:36:47Z raasch
42	! file re-formatted to follow the PALM coding standard
43	!
44	! 4560 2020-06-11 12:19:47Z suehring
45	! - Bugfix in calculation of vertical momentum and scalar fluxes
46	! - remove averaged output variables from PUBLIC list
47	!
48	! 4535 2020-05-15 12:07:23Z raasch
49	! bugfix for restart data format query
50	!
51	! 4518 2020-05-04 15:44:28Z suehring
52	! * Define arrays over ghost points in order to allow for standard mpi-io treatment. By this
53	! modularization of restart-data input is possible with the module interface.
54	! * Bugfix: add missing restart input of wtheta_av, wq_av, wu_av, and wv_av.
55	!
56	! 4517 2020-05-03 14:29:30Z raasch
57	! use statement for exchange horiz added,
58	! bugfix for call of exchange horiz 2d
59	!
60	! 4431 2020-02-27 23:23:01Z gronemeier
61	! added wspeed and wdir output; bugfix: set fill_value in case of masked output
62	!
63	! 4360 2020-01-07 11:25:50Z suehring
64	! added output of wu, wv, wtheta and wq to enable covariance calculation according to temporal EC
65	! method
66	!
67	! 4346 2019-12-18 11:55:56Z motisi
68	! Introduction of wall_flags_total_0, which currently sets bits based on static topography
69	! information used in wall_flags_static_0
70	!
71	! 4331 2019-12-10 18:25:02Z suehring
72	! - Modularize 2-m potential temperature output
73	! - New output for 10-m wind speed
74	!
75	! 4329 2019-12-10 15:46:36Z motisi
76	! Renamed wall_flags_0 to wall_flags_static_0
77	!
78	! 4182 2019-08-22 15:20:23Z scharf
79	! Corrected "Former revisions" section
80	!
81	! 4167 2019-08-16 11:01:48Z suehring
82	! Changed behaviour of masked output over surface to follow terrain and ignore buildings
83	! (J.Resler, T.Gronemeier)
84	!
85	! 4157 2019-08-14 09:19:12Z suehring
86	! Initialization restructured, in order to work also when data output during spin-up is enabled.
87	!
88	! 4132 2019-08-02 12:34:17Z suehring
89	! Bugfix in masked data output
90	!
91	! 4069 2019-07-01 14:05:51Z Giersch
92	! Masked output running index mid has been introduced as a local variable to avoid runtime error
93	! (Loop variable has been modified) in time_integration
94	!
95	! 4039 2019-06-18 10:32:41Z suehring
96	! - Add output of uu, vv, ww to enable variance calculation according temporal EC method
97	! - Allocate arrays only when they are required
98	! - Formatting adjustment
99	! - Rename subroutines
100	! - Further modularization
101	!
102	! 3998 2019-05-23 13:38:11Z suehring
103	! Bugfix in gathering all output strings
104	!
105	! 3995 2019-05-22 18:59:54Z suehring
106	! Avoid compiler warnings about unused variable and fix string operation which is not allowed with
107	! PGI compiler
108	!
109	! 3994 2019-05-22 18:08:09Z suehring
110	! Initial revision
111	!
112	! Authors:
113	! --------
114	! @author Farah Kanani-Suehring
115	!
116	!
117	! Description:
118	! ------------
119	!> ...
120	!--------------------------------------------------------------------------------------------------!
121	MODULE diagnostic_output_quantities_mod
122
123	USE arrays_3d, &
124	ONLY: ddzu, &
125	exner, &
126	hyp, &
127	pt, &
128	q, &
129	s, &
130	u, &
131	v, &
132	w, &
133	zu, &
134	zw
135
136	USE basic_constants_and_equations_mod, &
137	ONLY: degc_to_k, &
138	kappa, &
139	pi, &
140	magnus, &
141	rd_d_rv
142
143	USE control_parameters, &
144	ONLY: current_timestep_number, &
145	data_output, &
146	data_output_pr, &
147	humidity, &
148	message_string, &
149	passive_scalar, &
150	restart_data_format_output, &
151	varnamelength
152	!
153	! USE cpulog, &
154	! ONLY: cpu_log, log_point
155
156	USE exchange_horiz_mod, &
157	ONLY: exchange_horiz_2d
158
159	USE grid_variables, &
160	ONLY: ddx, ddy
161
162	USE indices, &
163	ONLY: nbgp, &
164	nxl, &
165	nxlg, &
166	nxr, &
167	nxrg, &
168	nyn, &
169	nyng, &
170	nys, &
171	nysg, &
172	nzb, &
173	nzt, &
174	wall_flags_total_0
175
176	USE kinds
177
178	USE profil_parameter, &
179	ONLY: dopr_index
180
181	USE restart_data_mpi_io_mod, &
182	ONLY: rd_mpi_io_check_array, &
183	rrd_mpi_io, &
184	wrd_mpi_io
185
186	USE statistics, &
187	ONLY: hom, &
188	pr_palm, &
189	rmask, &
190	statistic_regions, &
191	sums_l
192
193	USE surface_mod, &
194	ONLY: surf_def_h, &
195	surf_lsm_h, &
196	surf_type, &
197	surf_usm_h
198
199	IMPLICIT NONE
200
201	CHARACTER(LEN=varnamelength), DIMENSION(500) :: do_all = ' '
202
203	INTEGER(iwp) :: timestep_number_at_prev_calc = 0 !< ...at previous diagnostic output calculation
204
205	LOGICAL :: initialized_diagnostic_output_quantities = .FALSE. !< flag indicating whether output is initialized
206	LOGICAL :: prepared_diagnostic_output_quantities = .FALSE. !< flag indicating whether output is p
207
208	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pt_2m !< 2-m air potential temperature
209	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pt_2m_av !< averaged 2-m air potential temperature
210	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: uv_10m !< horizontal wind speed at 10m
211	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: uv_10m_av !< averaged horizontal wind speed at 10m
212
213	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: rh !< relative humidity
214	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: rh_av !< avg. relative humidity
215	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ta !< air temperature
216	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ta_av !< avg. air temperature
217	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ti !< rotation(u,v,w) aka turbulence intensity
218	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ti_av !< avg. rotation(u,v,w) aka turbulence intensity
219	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: u_center !< u at center of grid box
220	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: u_center_av !< mean of u_center
221	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: uu !< uu
222	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: uu_av !< mean of uu
223	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: wspeed !< horizontal wind speed
224	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: wspeed_av !< mean of horizotal wind speed
225	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: v_center !< v at center of grid box
226	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: v_center_av !< mean of v_center
227	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: vv !< vv
228	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: vv_av !< mean of vv
229	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: wdir !< wind direction
230	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: wdir_av !< mean wind direction
231	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ww !< ww
232	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ww_av !< mean of ww
233	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: wu !< wu
234	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: wu_av !< mean of wu
235	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: wv !< wv
236	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: wv_av !< mean of wv
237	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: wtheta !< wtheta
238	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: wtheta_av !< mean of wtheta
239	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: wq !< wq
240	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: wq_av !< mean of wq
241	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ws !< ws
242	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ws_av !< mean of ws
243
244
245	SAVE
246
247	PRIVATE
248
249	!
250	!-- Public variables
251	PUBLIC do_all, &
252	initialized_diagnostic_output_quantities, &
253	prepared_diagnostic_output_quantities, &
254	timestep_number_at_prev_calc
255	!
256	!-- Public routines
257	PUBLIC doq_3d_data_averaging, &
258	doq_calculate, &
259	doq_check_data_output, &
260	doq_check_data_output_pr, &
261	doq_define_netcdf_grid, &
262	doq_init, &
263	doq_output_2d, &
264	doq_output_3d, &
265	doq_output_mask, &
266	doq_statistics, &
267	doq_rrd_local, &
268	doq_wrd_local
269
270
271	INTERFACE doq_3d_data_averaging
272	MODULE PROCEDURE doq_3d_data_averaging
273	END INTERFACE doq_3d_data_averaging
274
275	INTERFACE doq_calculate
276	MODULE PROCEDURE doq_calculate
277	END INTERFACE doq_calculate
278
279	INTERFACE doq_check_data_output
280	MODULE PROCEDURE doq_check_data_output
281	END INTERFACE doq_check_data_output
282
283	INTERFACE doq_check_data_output_pr
284	MODULE PROCEDURE doq_check_data_output_pr
285	END INTERFACE doq_check_data_output_pr
286
287	INTERFACE doq_define_netcdf_grid
288	MODULE PROCEDURE doq_define_netcdf_grid
289	END INTERFACE doq_define_netcdf_grid
290
291	INTERFACE doq_output_2d
292	MODULE PROCEDURE doq_output_2d
293	END INTERFACE doq_output_2d
294
295	INTERFACE doq_output_3d
296	MODULE PROCEDURE doq_output_3d
297	END INTERFACE doq_output_3d
298
299	INTERFACE doq_output_mask
300	MODULE PROCEDURE doq_output_mask
301	END INTERFACE doq_output_mask
302
303	INTERFACE doq_init
304	MODULE PROCEDURE doq_init
305	END INTERFACE doq_init
306
307	INTERFACE doq_statistics
308	MODULE PROCEDURE doq_statistics
309	END INTERFACE doq_statistics
310
311	INTERFACE doq_prepare
312	MODULE PROCEDURE doq_prepare
313	END INTERFACE doq_prepare
314
315	INTERFACE doq_rrd_local
316	MODULE PROCEDURE doq_rrd_local_ftn
317	MODULE PROCEDURE doq_rrd_local_mpi
318	END INTERFACE doq_rrd_local
319
320	INTERFACE doq_wrd_local
321	MODULE PROCEDURE doq_wrd_local
322	END INTERFACE doq_wrd_local
323
324	CONTAINS
325
326	!--------------------------------------------------------------------------------------------------!
327	! Description:
328	! ------------
329	!> Sum up and time-average diagnostic output quantities as well as allocate the array necessary for
330	!> storing the average.
331	!--------------------------------------------------------------------------------------------------!
332	SUBROUTINE doq_3d_data_averaging( mode, variable )
333
334	USE control_parameters, &
335	ONLY: average_count_3d
336
337	CHARACTER (LEN=*) :: mode !<
338	CHARACTER (LEN=*) :: variable !<
339
340	INTEGER(iwp) :: i !<
341	INTEGER(iwp) :: j !<
342	INTEGER(iwp) :: k !<
343
344
345	IF ( mode == 'allocate' ) THEN
346
347	SELECT CASE ( TRIM( variable ) )
348
349	CASE ( 'rh' )
350	IF ( .NOT. ALLOCATED( rh_av ) ) THEN
351	ALLOCATE( rh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
352	ENDIF
353	rh_av = 0.0_wp
354
355	CASE ( 'ta' )
356	IF ( .NOT. ALLOCATED( ta_av ) ) THEN
357	ALLOCATE( ta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
358	ENDIF
359	ta_av = 0.0_wp
360
361	CASE ( 'ti' )
362	IF ( .NOT. ALLOCATED( ti_av ) ) THEN
363	ALLOCATE( ti_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
364	ENDIF
365	ti_av = 0.0_wp
366
367	CASE ( 'uu' )
368	IF ( .NOT. ALLOCATED( uu_av ) ) THEN
369	ALLOCATE( uu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
370	ENDIF
371	uu_av = 0.0_wp
372
373	CASE ( 'vv' )
374	IF ( .NOT. ALLOCATED( vv_av ) ) THEN
375	ALLOCATE( vv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
376	ENDIF
377	vv_av = 0.0_wp
378
379	CASE ( 'ww' )
380	IF ( .NOT. ALLOCATED( ww_av ) ) THEN
381	ALLOCATE( ww_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
382	ENDIF
383	ww_av = 0.0_wp
384
385	CASE ( 'wu' )
386	IF ( .NOT. ALLOCATED( wu_av ) ) THEN
387	ALLOCATE( wu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
388	ENDIF
389	wu_av = 0.0_wp
390
391	CASE ( 'wv' )
392	IF ( .NOT. ALLOCATED( wv_av ) ) THEN
393	ALLOCATE( wv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
394	ENDIF
395	wv_av = 0.0_wp
396
397	CASE ( 'wtheta' )
398	IF ( .NOT. ALLOCATED( wtheta_av ) ) THEN
399	ALLOCATE( wtheta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
400	ENDIF
401	wtheta_av = 0.0_wp
402
403	CASE ( 'wq' )
404	IF ( .NOT. ALLOCATED( wq_av ) ) THEN
405	ALLOCATE( wq_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
406	ENDIF
407	wq_av = 0.0_wp
408
409	CASE ( 'ws' )
410	IF ( .NOT. ALLOCATED( ws_av ) ) THEN
411	ALLOCATE( ws_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
412	ENDIF
413	ws_av = 0.0_wp
414
415	CASE ( 'theta_2m*' )
416	IF ( .NOT. ALLOCATED( pt_2m_av ) ) THEN
417	ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )
418	ENDIF
419	pt_2m_av = 0.0_wp
420
421	CASE ( 'wspeed_10m*' )
422	IF ( .NOT. ALLOCATED( uv_10m_av ) ) THEN
423	ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) )
424	ENDIF
425	uv_10m_av = 0.0_wp
426
427	CASE ( 'wspeed' )
428	IF ( .NOT. ALLOCATED( wspeed_av ) ) THEN
429	ALLOCATE( wspeed_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
430	ENDIF
431	wspeed_av = 0.0_wp
432
433	CASE ( 'wdir' )
434	IF ( .NOT. ALLOCATED( u_center_av ) ) THEN
435	ALLOCATE( u_center_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
436	ENDIF
437	IF ( .NOT. ALLOCATED( v_center_av ) ) THEN
438	ALLOCATE( v_center_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
439	ENDIF
440	u_center_av = 0.0_wp
441	v_center_av = 0.0_wp
442
443	CASE DEFAULT
444	CONTINUE
445
446	END SELECT
447
448	ELSEIF ( mode == 'sum' ) THEN
449
450	SELECT CASE ( TRIM( variable ) )
451
452	CASE ( 'rh' )
453	IF ( ALLOCATED( rh_av ) ) THEN
454	DO i = nxl, nxr
455	DO j = nys, nyn
456	DO k = nzb, nzt+1
457	rh_av(k,j,i) = rh_av(k,j,i) + rh(k,j,i)
458	ENDDO
459	ENDDO
460	ENDDO
461	ENDIF
462
463	CASE ( 'ta' )
464	IF ( ALLOCATED( ta_av ) ) THEN
465	DO i = nxl, nxr
466	DO j = nys, nyn
467	DO k = nzb, nzt+1
468	ta_av(k,j,i) = ta_av(k,j,i) + ta(k,j,i)
469	ENDDO
470	ENDDO
471	ENDDO
472	ENDIF
473
474	CASE ( 'ti' )
475	IF ( ALLOCATED( ti_av ) ) THEN
476	DO i = nxl, nxr
477	DO j = nys, nyn
478	DO k = nzb, nzt+1
479	ti_av(k,j,i) = ti_av(k,j,i) + ti(k,j,i)
480	ENDDO
481	ENDDO
482	ENDDO
483	ENDIF
484
485	CASE ( 'uu' )
486	IF ( ALLOCATED( uu_av ) ) THEN
487	DO i = nxl, nxr
488	DO j = nys, nyn
489	DO k = nzb, nzt+1
490	uu_av(k,j,i) = uu_av(k,j,i) + uu(k,j,i)
491	ENDDO
492	ENDDO
493	ENDDO
494	ENDIF
495
496	CASE ( 'vv' )
497	IF ( ALLOCATED( vv_av ) ) THEN
498	DO i = nxl, nxr
499	DO j = nys, nyn
500	DO k = nzb, nzt+1
501	vv_av(k,j,i) = vv_av(k,j,i) + vv(k,j,i)
502	ENDDO
503	ENDDO
504	ENDDO
505	ENDIF
506
507	CASE ( 'ww' )
508	IF ( ALLOCATED( ww_av ) ) THEN
509	DO i = nxl, nxr
510	DO j = nys, nyn
511	DO k = nzb, nzt+1
512	ww_av(k,j,i) = ww_av(k,j,i) + ww(k,j,i)
513	ENDDO
514	ENDDO
515	ENDDO
516	ENDIF
517
518	CASE ( 'wu' )
519	IF ( ALLOCATED( wu_av ) ) THEN
520	DO i = nxl, nxr
521	DO j = nys, nyn
522	DO k = nzb, nzt+1
523	wu_av(k,j,i) = wu_av(k,j,i) + wu(k,j,i)
524	ENDDO
525	ENDDO
526	ENDDO
527	ENDIF
528
529	CASE ( 'wv' )
530	IF ( ALLOCATED( wv_av ) ) THEN
531	DO i = nxl, nxr
532	DO j = nys, nyn
533	DO k = nzb, nzt+1
534	wv_av(k,j,i) = wv_av(k,j,i) + wv(k,j,i)
535	ENDDO
536	ENDDO
537	ENDDO
538	ENDIF
539
540	CASE ( 'wtheta' )
541	IF ( ALLOCATED( wtheta_av ) ) THEN
542	DO i = nxl, nxr
543	DO j = nys, nyn
544	DO k = nzb, nzt+1
545	wtheta_av(k,j,i) = wtheta_av(k,j,i) + wtheta(k,j,i)
546	ENDDO
547	ENDDO
548	ENDDO
549	ENDIF
550
551	CASE ( 'wq' )
552	IF ( ALLOCATED( wq_av ) ) THEN
553	DO i = nxl, nxr
554	DO j = nys, nyn
555	DO k = nzb, nzt+1
556	wq_av(k,j,i) = wq_av(k,j,i) + wq(k,j,i)
557	ENDDO
558	ENDDO
559	ENDDO
560	ENDIF
561
562	CASE ( 'ws' )
563	IF ( ALLOCATED( ws_av ) ) THEN
564	DO i = nxl, nxr
565	DO j = nys, nyn
566	DO k = nzb, nzt+1
567	ws_av(k,j,i) = ws_av(k,j,i) + ws(k,j,i)
568	ENDDO
569	ENDDO
570	ENDDO
571	ENDIF
572
573	CASE ( 'theta_2m*' )
574	IF ( ALLOCATED( pt_2m_av ) ) THEN
575	DO i = nxl, nxr
576	DO j = nys, nyn
577	pt_2m_av(j,i) = pt_2m_av(j,i) + pt_2m(j,i)
578	ENDDO
579	ENDDO
580	ENDIF
581
582	CASE ( 'wspeed_10m*' )
583	IF ( ALLOCATED( uv_10m_av ) ) THEN
584	DO i = nxl, nxr
585	DO j = nys, nyn
586	uv_10m_av(j,i) = uv_10m_av(j,i) + uv_10m(j,i)
587	ENDDO
588	ENDDO
589	ENDIF
590
591	CASE ( 'wspeed' )
592	IF ( ALLOCATED( wspeed_av ) ) THEN
593	DO i = nxl, nxr
594	DO j = nys, nyn
595	DO k = nzb, nzt+1
596	wspeed_av(k,j,i) = wspeed_av(k,j,i) + wspeed(k,j,i)
597	ENDDO
598	ENDDO
599	ENDDO
600	ENDIF
601
602	CASE ( 'wdir' )
603	IF ( ALLOCATED( u_center_av ) .AND. ALLOCATED( v_center_av ) ) THEN
604	DO i = nxl, nxr
605	DO j = nys, nyn
606	DO k = nzb, nzt+1
607	u_center_av(k,j,i) = u_center_av(k,j,i) + u_center(k,j,i)
608	v_center_av(k,j,i) = v_center_av(k,j,i) + v_center(k,j,i)
609	ENDDO
610	ENDDO
611	ENDDO
612	ENDIF
613
614	CASE DEFAULT
615	CONTINUE
616
617	END SELECT
618
619	ELSEIF ( mode == 'average' ) THEN
620
621	SELECT CASE ( TRIM( variable ) )
622
623	CASE ( 'rh' )
624	IF ( ALLOCATED( rh_av ) ) THEN
625	DO i = nxl, nxr
626	DO j = nys, nyn
627	DO k = nzb, nzt+1
628	rh_av(k,j,i) = rh_av(k,j,i) / REAL( average_count_3d, KIND=wp )
629	ENDDO
630	ENDDO
631	ENDDO
632	ENDIF
633
634	CASE ( 'ta' )
635	IF ( ALLOCATED( ta_av ) ) THEN
636	DO i = nxl, nxr
637	DO j = nys, nyn
638	DO k = nzb, nzt+1
639	ta_av(k,j,i) = ta_av(k,j,i) / REAL( average_count_3d, KIND=wp )
640	ENDDO
641	ENDDO
642	ENDDO
643	ENDIF
644
645	CASE ( 'ti' )
646	IF ( ALLOCATED( ti_av ) ) THEN
647	DO i = nxl, nxr
648	DO j = nys, nyn
649	DO k = nzb, nzt+1
650	ti_av(k,j,i) = ti_av(k,j,i) / REAL( average_count_3d, KIND=wp )
651	ENDDO
652	ENDDO
653	ENDDO
654	ENDIF
655
656	CASE ( 'uu' )
657	IF ( ALLOCATED( uu_av ) ) THEN
658	DO i = nxl, nxr
659	DO j = nys, nyn
660	DO k = nzb, nzt+1
661	uu_av(k,j,i) = uu_av(k,j,i) / REAL( average_count_3d, KIND=wp )
662	ENDDO
663	ENDDO
664	ENDDO
665	ENDIF
666
667	CASE ( 'vv' )
668	IF ( ALLOCATED( vv_av ) ) THEN
669	DO i = nxl, nxr
670	DO j = nys, nyn
671	DO k = nzb, nzt+1
672	vv_av(k,j,i) = vv_av(k,j,i) / REAL( average_count_3d, KIND=wp )
673	ENDDO
674	ENDDO
675	ENDDO
676	ENDIF
677
678	CASE ( 'ww' )
679	IF ( ALLOCATED( ww_av ) ) THEN
680	DO i = nxl, nxr
681	DO j = nys, nyn
682	DO k = nzb, nzt+1
683	ww_av(k,j,i) = ww_av(k,j,i) / REAL( average_count_3d, KIND=wp )
684	ENDDO
685	ENDDO
686	ENDDO
687	ENDIF
688
689	CASE ( 'wu' )
690	IF ( ALLOCATED( wu_av ) ) THEN
691	DO i = nxl, nxr
692	DO j = nys, nyn
693	DO k = nzb, nzt+1
694	wu_av(k,j,i) = wu_av(k,j,i) / REAL( average_count_3d, KIND=wp )
695	ENDDO
696	ENDDO
697	ENDDO
698	ENDIF
699
700	CASE ( 'wv' )
701	IF ( ALLOCATED( wv_av ) ) THEN
702	DO i = nxl, nxr
703	DO j = nys, nyn
704	DO k = nzb, nzt+1
705	wv_av(k,j,i) = wv_av(k,j,i) / REAL( average_count_3d, KIND=wp )
706	ENDDO
707	ENDDO
708	ENDDO
709	ENDIF
710
711	CASE ( 'wtheta' )
712	IF ( ALLOCATED( wtheta_av ) ) THEN
713	DO i = nxl, nxr
714	DO j = nys, nyn
715	DO k = nzb, nzt+1
716	wtheta_av(k,j,i) = wtheta_av(k,j,i) / REAL( average_count_3d, KIND=wp )
717	ENDDO
718	ENDDO
719	ENDDO
720	ENDIF
721
722	CASE ( 'wq' )
723	IF ( ALLOCATED( wq_av ) ) THEN
724	DO i = nxl, nxr
725	DO j = nys, nyn
726	DO k = nzb, nzt+1
727	wq_av(k,j,i) = wq_av(k,j,i) / REAL( average_count_3d, KIND=wp )
728	ENDDO
729	ENDDO
730	ENDDO
731	ENDIF
732
733	CASE ( 'ws' )
734	IF ( ALLOCATED( ws_av ) ) THEN
735	DO i = nxl, nxr
736	DO j = nys, nyn
737	DO k = nzb, nzt+1
738	ws_av(k,j,i) = ws_av(k,j,i) / REAL( average_count_3d, KIND=wp )
739	ENDDO
740	ENDDO
741	ENDDO
742	ENDIF
743
744	CASE ( 'theta_2m*' )
745	IF ( ALLOCATED( pt_2m_av ) ) THEN
746	DO i = nxlg, nxrg
747	DO j = nysg, nyng
748	pt_2m_av(j,i) = pt_2m_av(j,i) / REAL( average_count_3d, KIND=wp )
749	ENDDO
750	ENDDO
751	CALL exchange_horiz_2d( pt_2m_av )
752	ENDIF
753
754	CASE ( 'wspeed_10m*' )
755	IF ( ALLOCATED( uv_10m_av ) ) THEN
756	DO i = nxlg, nxrg
757	DO j = nysg, nyng
758	uv_10m_av(j,i) = uv_10m_av(j,i) / REAL( average_count_3d, KIND=wp )
759	ENDDO
760	ENDDO
761	CALL exchange_horiz_2d( uv_10m_av )
762	ENDIF
763
764	CASE ( 'wspeed' )
765	IF ( ALLOCATED( wspeed_av ) ) THEN
766	DO i = nxl, nxr
767	DO j = nys, nyn
768	DO k = nzb, nzt+1
769	wspeed_av(k,j,i) = wspeed_av(k,j,i) / REAL( average_count_3d, KIND=wp )
770	ENDDO
771	ENDDO
772	ENDDO
773	ENDIF
774
775	CASE ( 'wdir' )
776	IF ( ALLOCATED( u_center_av ) .AND. ALLOCATED( v_center_av ) ) THEN
777
778	IF ( .NOT. ALLOCATED( wdir_av ) ) THEN
779	ALLOCATE( wdir_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
780	ENDIF
781	wdir_av = 0.0_wp
782
783	DO i = nxl, nxr
784	DO j = nys, nyn
785	DO k = nzb, nzt+1
786	u_center_av(k,j,i) = u_center_av(k,j,i) / REAL( average_count_3d, KIND=wp )
787	v_center_av(k,j,i) = v_center_av(k,j,i) / REAL( average_count_3d, KIND=wp )
788	wdir_av(k,j,i) = ATAN2( u_center_av(k,j,i), v_center_av(k,j,i) ) &
789	/ pi * 180.0_wp + 180.0_wp
790	ENDDO
791	ENDDO
792	ENDDO
793	ENDIF
794
795	END SELECT
796
797	ENDIF
798
799
800	END SUBROUTINE doq_3d_data_averaging
801
802	!--------------------------------------------------------------------------------------------------!
803	! Description:
804	! ------------
805	!> Check data output for diagnostic output
806	!--------------------------------------------------------------------------------------------------!
807	SUBROUTINE doq_check_data_output( var, unit, i, ilen, k )
808
809	IMPLICIT NONE
810
811	CHARACTER (LEN=*) :: unit !<
812	CHARACTER (LEN=*) :: var !<
813
814	INTEGER(iwp), OPTIONAL, INTENT(IN) :: i !< Current element of data_output
815	INTEGER(iwp), OPTIONAL, INTENT(IN) :: ilen !< Length of current entry in data_output
816	INTEGER(iwp), OPTIONAL, INTENT(IN) :: k !< Output is xy mode? 0 = no, 1 = yes
817
818
819	SELECT CASE ( TRIM( var ) )
820
821	CASE ( 'rh' )
822	unit = '%'
823
824	CASE ( 'ta' )
825	unit = 'degree_C'
826
827	CASE ( 'ti' )
828	unit = '1/s'
829
830	CASE ( 'uu' )
831	unit = 'm2/s2'
832
833	CASE ( 'vv' )
834	unit = 'm2/s2'
835
836	CASE ( 'ww' )
837	unit = 'm2/s2'
838
839	CASE ( 'wu' )
840	unit = 'm2/s2'
841
842	CASE ( 'wv' )
843	unit = 'm2/s2'
844
845	CASE ( 'wtheta' )
846	unit = 'Km/s'
847
848	CASE ( 'wq' )
849	unit = 'm/s'
850
851	CASE ( 'ws' )
852	IF ( .NOT. passive_scalar ) THEN
853	message_string = 'data_output = ' // TRIM( var ) // &
854	' is not implemented for passive_scalar = .FALSE.'
855	CALL message( 'diagnostic_output', 'PA0093', 1, 2, 0, 6, 0 )
856	ENDIF
857	unit = 'm/s'
858
859	CASE ( 'wspeed' )
860	unit = 'm/s'
861
862	CASE ( 'wdir' )
863	unit = 'degree'
864	!
865	!-- Treat horizotal cross-section output quanatities
866	CASE ( 'theta_2m', 'wspeed_10m' )
867	!
868	!-- Check if output quantity is _xy only.
869	IF ( k == 0 .OR. data_output(i)(ilen-2:ilen) /= '_xy' ) THEN
870	message_string = 'illegal value for data_output: "' // &
871	TRIM( var ) // '" & only 2d-horizontal ' // &
872	'cross sections are allowed for this value'
873	CALL message( 'diagnostic_output', 'PA0111', 1, 2, 0, 6, 0 )
874	ENDIF
875
876	IF ( TRIM( var ) == 'theta_2m*' ) unit = 'K'
877	IF ( TRIM( var ) == 'wspeed_10m*' ) unit = 'm/s'
878
879	CASE DEFAULT
880	unit = 'illegal'
881
882	END SELECT
883
884
885	END SUBROUTINE doq_check_data_output
886
887
888	!--------------------------------------------------------------------------------------------------!
889	! Description:
890	! ------------
891	!> Set the unit of user defined profile output quantities. For those variables not recognized by the
892	!> user, the parameter unit is set to "illegal", which tells the calling routine that the
893	!> output variable is not defined and leads to a program abort.
894	!--------------------------------------------------------------------------------------------------!
895	SUBROUTINE doq_check_data_output_pr( variable, var_count, unit, dopr_unit )
896
897	CHARACTER (LEN=*) :: unit !<
898	CHARACTER (LEN=*) :: variable !<
899	CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit
900
901	INTEGER(iwp) :: pr_index !<
902	INTEGER(iwp) :: var_count !<
903
904	SELECT CASE ( TRIM( variable ) )
905
906	CASE ( 'rh' )
907	IF ( .NOT. humidity ) THEN
908	message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) // &
909	' requires humidity'
910	CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
911	ENDIF
912	pr_index = 130
913	dopr_index(var_count) = pr_index
914	dopr_unit = '%'
915	unit = dopr_unit
916	hom(:,2,pr_index,:) = SPREAD( zu, 2, statistic_regions+1 )
917
918	CASE DEFAULT
919	unit = 'illegal'
920
921	END SELECT
922
923
924	END SUBROUTINE doq_check_data_output_pr
925
926
927	!--------------------------------------------------------------------------------------------------!
928	!
929	! Description:
930	! ------------
931	!> Subroutine defining appropriate grid for netcdf variables.
932	!--------------------------------------------------------------------------------------------------!
933	SUBROUTINE doq_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z )
934
935	IMPLICIT NONE
936
937	CHARACTER (LEN=*), INTENT(IN) :: variable !<
938	CHARACTER (LEN=*), INTENT(OUT) :: grid_x !<
939	CHARACTER (LEN=*), INTENT(OUT) :: grid_y !<
940	CHARACTER (LEN=*), INTENT(OUT) :: grid_z !<
941
942	LOGICAL, INTENT(OUT) :: found !<
943
944
945	found = .TRUE.
946
947	SELECT CASE ( TRIM( variable ) )
948	!
949	!-- s grid
950	CASE ( 'rh', 'rh_xy', 'rh_xz', 'rh_yz', &
951	'ta', 'ta_xy', 'ta_xz', 'ta_yz', &
952	'ti', 'ti_xy', 'ti_xz', 'ti_yz', &
953	'wspeed', 'wspeed_xy', 'wspeed_xz', 'wspeed_yz', &
954	'wdir', 'wdir_xy', 'wdir_xz', 'wdir_yz' )
955
956	grid_x = 'x'
957	grid_y = 'y'
958	grid_z = 'zu'
959	!
960	!-- s grid surface variables
961	CASE ( 'theta_2m_xy', 'wspeed_10m' )
962
963	grid_x = 'x'
964	grid_y = 'y'
965	grid_z = 'zu'
966	!
967	!-- u grid
968	CASE ( 'uu', 'uu_xy', 'uu_xz', 'uu_yz' )
969
970	grid_x = 'xu'
971	grid_y = 'y'
972	grid_z = 'zu'
973	!
974	!-- v grid
975	CASE ( 'vv', 'vv_xy', 'vv_xz', 'vv_yz' )
976
977	grid_x = 'x'
978	grid_y = 'yv'
979	grid_z = 'zu'
980	!
981	!-- w grid
982	CASE ( 'ww', 'ww_xy', 'ww_xz', 'ww_yz', &
983	'wu', 'wu_xy', 'wu_xz', 'wu_yz', &
984	'wv', 'wv_xy', 'wv_xz', 'wv_yz', &
985	'wtheta', 'wtheta_xy', 'wtheta_xz', 'wtheta_yz', &
986	'wq', 'wq_xy', 'wq_xz', 'wq_yz', &
987	'ws', 'ws_xy', 'ws_xz', 'ws_yz' )
988
989	grid_x = 'x'
990	grid_y = 'y'
991	grid_z = 'zw'
992
993	CASE DEFAULT
994	found = .FALSE.
995	grid_x = 'none'
996	grid_y = 'none'
997	grid_z = 'none'
998
999	END SELECT
1000
1001
1002	END SUBROUTINE doq_define_netcdf_grid
1003
1004	!--------------------------------------------------------------------------------------------------!
1005	!
1006	! Description:
1007	! ------------
1008	!> Subroutine defining 2D output variables
1009	!--------------------------------------------------------------------------------------------------!
1010	SUBROUTINE doq_output_2d( av, variable, found, grid, mode, local_pf, two_d, nzb_do, nzt_do, &
1011	fill_value )
1012
1013
1014	IMPLICIT NONE
1015
1016	CHARACTER (LEN=*) :: grid !<
1017	CHARACTER (LEN=*) :: mode !<
1018	CHARACTER (LEN=*) :: variable !<
1019
1020	INTEGER(iwp) :: av !< value indicating averaged or non-averaged output
1021	INTEGER(iwp) :: flag_nr !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
1022	INTEGER(iwp) :: i !< grid index x-direction
1023	INTEGER(iwp) :: j !< grid index y-direction
1024	INTEGER(iwp) :: k !< grid index z-direction
1025	INTEGER(iwp) :: nzb_do !<
1026	INTEGER(iwp) :: nzt_do !<
1027
1028	LOGICAL :: found !< true if variable is in list
1029	LOGICAL :: resorted !< true if array is resorted
1030	LOGICAL :: two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
1031
1032	REAL(wp) :: fill_value !< value for the _FillValue attribute
1033
1034	REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !<
1035	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which needs to be resorted for output
1036
1037
1038	flag_nr = 0
1039	found = .TRUE.
1040	resorted = .FALSE.
1041	two_d = .FALSE.
1042
1043	SELECT CASE ( TRIM( variable ) )
1044
1045	CASE ( 'rh_xy', 'rh_xz', 'rh_yz' )
1046	IF ( av == 0 ) THEN
1047	to_be_resorted => rh
1048	ELSE
1049	IF ( .NOT. ALLOCATED( rh_av ) ) THEN
1050	ALLOCATE( rh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1051	rh_av = REAL( fill_value, KIND = wp )
1052	ENDIF
1053	to_be_resorted => rh_av
1054	ENDIF
1055	flag_nr = 0
1056
1057	IF ( mode == 'xy' ) grid = 'zu'
1058
1059	CASE ( 'ta_xy', 'ta_xz', 'ta_yz' )
1060	IF ( av == 0 ) THEN
1061	to_be_resorted => ta
1062	ELSE
1063	IF ( .NOT. ALLOCATED( ta_av ) ) THEN
1064	ALLOCATE( ta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1065	ta_av = REAL( fill_value, KIND = wp )
1066	ENDIF
1067	to_be_resorted => ta_av
1068	ENDIF
1069	flag_nr = 0
1070
1071	IF ( mode == 'xy' ) grid = 'zu'
1072
1073	CASE ( 'ti_xy', 'ti_xz', 'ti_yz' )
1074	IF ( av == 0 ) THEN
1075	to_be_resorted => ti
1076	ELSE
1077	IF ( .NOT. ALLOCATED( ti_av ) ) THEN
1078	ALLOCATE( ti_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1079	ti_av = REAL( fill_value, KIND = wp )
1080	ENDIF
1081	to_be_resorted => ti_av
1082	ENDIF
1083	flag_nr = 0
1084
1085	IF ( mode == 'xy' ) grid = 'zu'
1086
1087	CASE ( 'uu_xy', 'uu_xz', 'uu_yz' )
1088	IF ( av == 0 ) THEN
1089	to_be_resorted => uu
1090	ELSE
1091	IF ( .NOT. ALLOCATED( uu_av ) ) THEN
1092	ALLOCATE( uu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1093	uu_av = REAL( fill_value, KIND = wp )
1094	ENDIF
1095	to_be_resorted => uu_av
1096	ENDIF
1097	flag_nr = 1
1098
1099	IF ( mode == 'xy' ) grid = 'zu'
1100
1101	CASE ( 'vv_xy', 'vv_xz', 'vv_yz' )
1102	IF ( av == 0 ) THEN
1103	to_be_resorted => vv
1104	ELSE
1105	IF ( .NOT. ALLOCATED( vv_av ) ) THEN
1106	ALLOCATE( vv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1107	vv_av = REAL( fill_value, KIND = wp )
1108	ENDIF
1109	to_be_resorted => vv_av
1110	ENDIF
1111	flag_nr = 2
1112
1113	IF ( mode == 'xy' ) grid = 'zu'
1114
1115	CASE ( 'ww_xy', 'ww_xz', 'ww_yz' )
1116	IF ( av == 0 ) THEN
1117	to_be_resorted => ww
1118	ELSE
1119	IF ( .NOT. ALLOCATED( ww_av ) ) THEN
1120	ALLOCATE( ww_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1121	ww_av = REAL( fill_value, KIND = wp )
1122	ENDIF
1123	to_be_resorted => ww_av
1124	ENDIF
1125	flag_nr = 3
1126
1127	IF ( mode == 'xy' ) grid = 'zw'
1128
1129	CASE ( 'wu_xy', 'wu_xz', 'wu_yz' )
1130	IF ( av == 0 ) THEN
1131	to_be_resorted => wu
1132	ELSE
1133	IF ( .NOT. ALLOCATED( wu_av ) ) THEN
1134	ALLOCATE( wu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1135	wu_av = REAL( fill_value, KIND = wp )
1136	ENDIF
1137	to_be_resorted => wu_av
1138	ENDIF
1139	flag_nr = 0
1140
1141	IF ( mode == 'xy' ) grid = 'zw'
1142
1143	CASE ( 'wv_xy', 'wv_xz', 'wv_yz' )
1144	IF ( av == 0 ) THEN
1145	to_be_resorted => wv
1146	ELSE
1147	IF ( .NOT. ALLOCATED( wv_av ) ) THEN
1148	ALLOCATE( wv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1149	wv_av = REAL( fill_value, KIND = wp )
1150	ENDIF
1151	to_be_resorted => wv_av
1152	ENDIF
1153	flag_nr = 0
1154
1155	IF ( mode == 'xy' ) grid = 'zw'
1156
1157	CASE ( 'wtheta_xy', 'wtheta_xz', 'wtheta_yz' )
1158	IF ( av == 0 ) THEN
1159	to_be_resorted => wtheta
1160	ELSE
1161	IF ( .NOT. ALLOCATED( wtheta_av ) ) THEN
1162	ALLOCATE( wtheta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1163	wtheta_av = REAL( fill_value, KIND = wp )
1164	ENDIF
1165	to_be_resorted => wtheta_av
1166	ENDIF
1167	flag_nr = 0
1168
1169	IF ( mode == 'xy' ) grid = 'zw'
1170
1171	CASE ( 'wq_xy', 'wq_xz', 'wq_yz' )
1172	IF ( av == 0 ) THEN
1173	to_be_resorted => wq
1174	ELSE
1175	IF ( .NOT. ALLOCATED( wq_av ) ) THEN
1176	ALLOCATE( wq_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1177	wq_av = REAL( fill_value, KIND = wp )
1178	ENDIF
1179	to_be_resorted => wq_av
1180	ENDIF
1181	flag_nr = 0
1182
1183	IF ( mode == 'xy' ) grid = 'zw'
1184
1185	CASE ( 'ws_xy', 'ws_xz', 'ws_yz' )
1186	IF ( av == 0 ) THEN
1187	to_be_resorted => ws
1188	ELSE
1189	IF ( .NOT. ALLOCATED( ws_av ) ) THEN
1190	ALLOCATE( ws_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1191	ws_av = REAL( fill_value, KIND = wp )
1192	ENDIF
1193	to_be_resorted => ws_av
1194	ENDIF
1195	flag_nr = 0
1196
1197	IF ( mode == 'xy' ) grid = 'zw'
1198
1199	CASE ( 'theta_2m*_xy' ) ! 2d-array
1200	IF ( av == 0 ) THEN
1201	DO i = nxl, nxr
1202	DO j = nys, nyn
1203	local_pf(i,j,nzb+1) = pt_2m(j,i)
1204	ENDDO
1205	ENDDO
1206	ELSE
1207	IF ( .NOT. ALLOCATED( pt_2m_av ) ) THEN
1208	ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )
1209	pt_2m_av = REAL( fill_value, KIND = wp )
1210	ENDIF
1211	DO i = nxl, nxr
1212	DO j = nys, nyn
1213	local_pf(i,j,nzb+1) = pt_2m_av(j,i)
1214	ENDDO
1215	ENDDO
1216	ENDIF
1217	resorted = .TRUE.
1218	two_d = .TRUE.
1219	grid = 'zu1'
1220
1221	CASE ( 'wspeed_10m*_xy' ) ! 2d-array
1222	IF ( av == 0 ) THEN
1223	DO i = nxl, nxr
1224	DO j = nys, nyn
1225	local_pf(i,j,nzb+1) = uv_10m(j,i)
1226	ENDDO
1227	ENDDO
1228	ELSE
1229	IF ( .NOT. ALLOCATED( uv_10m_av ) ) THEN
1230	ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) )
1231	uv_10m_av = REAL( fill_value, KIND = wp )
1232	ENDIF
1233	DO i = nxl, nxr
1234	DO j = nys, nyn
1235	local_pf(i,j,nzb+1) = uv_10m_av(j,i)
1236	ENDDO
1237	ENDDO
1238	ENDIF
1239	resorted = .TRUE.
1240	two_d = .TRUE.
1241	grid = 'zu1'
1242
1243	CASE ( 'wspeed_xy', 'wspeed_xz', 'wspeed_yz' )
1244	IF ( av == 0 ) THEN
1245	to_be_resorted => wspeed
1246	ELSE
1247	IF ( .NOT. ALLOCATED( wspeed_av ) ) THEN
1248	ALLOCATE( wspeed_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1249	wspeed_av = REAL( fill_value, KIND = wp )
1250	ENDIF
1251	to_be_resorted => wspeed_av
1252	ENDIF
1253	flag_nr = 0
1254
1255	IF ( mode == 'xy' ) grid = 'zu'
1256
1257	CASE ( 'wdir_xy', 'wdir_xz', 'wdir_yz' )
1258	IF ( av == 0 ) THEN
1259	to_be_resorted => wdir
1260	ELSE
1261	IF ( .NOT. ALLOCATED( wdir_av ) ) THEN
1262	ALLOCATE( wdir_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1263	wdir_av = REAL( fill_value, KIND = wp )
1264	ENDIF
1265	to_be_resorted => wdir_av
1266	ENDIF
1267	flag_nr = 0
1268
1269	IF ( mode == 'xy' ) grid = 'zu'
1270
1271	CASE DEFAULT
1272	found = .FALSE.
1273	grid = 'none'
1274
1275	END SELECT
1276
1277	IF ( found .AND. .NOT. resorted ) THEN
1278	DO i = nxl, nxr
1279	DO j = nys, nyn
1280	DO k = nzb_do, nzt_do
1281	local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), &
1282	REAL( fill_value, KIND = wp ), &
1283	BTEST( wall_flags_total_0(k,j,i), flag_nr ) )
1284	ENDDO
1285	ENDDO
1286	ENDDO
1287	ENDIF
1288
1289	END SUBROUTINE doq_output_2d
1290
1291
1292	!--------------------------------------------------------------------------------------------------!
1293	!
1294	! Description:
1295	! ------------
1296	!> Subroutine defining 3D output variables
1297	!--------------------------------------------------------------------------------------------------!
1298	SUBROUTINE doq_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
1299
1300	IMPLICIT NONE
1301
1302	CHARACTER (LEN=*) :: variable !<
1303
1304	INTEGER(iwp) :: av !< index indicating averaged or instantaneous output
1305	INTEGER(iwp) :: flag_nr !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
1306	INTEGER(iwp) :: i !< index variable along x-direction
1307	INTEGER(iwp) :: j !< index variable along y-direction
1308	INTEGER(iwp) :: k !< index variable along z-direction
1309	INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0)
1310	INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d)
1311
1312	LOGICAL :: found !< true if variable is in list
1313	LOGICAL :: resorted !< true if array is resorted
1314
1315	REAL(wp) :: fill_value !< value for the _FillValue attribute
1316
1317	REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !<
1318	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which needs to be resorted for output
1319
1320
1321	flag_nr = 0
1322	found = .TRUE.
1323	resorted = .FALSE.
1324
1325	SELECT CASE ( TRIM( variable ) )
1326
1327	CASE ( 'rh' )
1328	IF ( av == 0 ) THEN
1329	to_be_resorted => rh
1330	ELSE
1331	IF ( .NOT. ALLOCATED( rh ) ) THEN
1332	ALLOCATE( rh(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1333	rh = REAL( fill_value, KIND = wp )
1334	ENDIF
1335	to_be_resorted => rh
1336	ENDIF
1337	flag_nr = 0
1338
1339	CASE ( 'ta' )
1340	IF ( av == 0 ) THEN
1341	to_be_resorted => ta
1342	ELSE
1343	IF ( .NOT. ALLOCATED( ta ) ) THEN
1344	ALLOCATE( ta(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1345	ta = REAL( fill_value, KIND = wp )
1346	ENDIF
1347	to_be_resorted => ta
1348	ENDIF
1349	flag_nr = 0
1350
1351	CASE ( 'ti' )
1352	IF ( av == 0 ) THEN
1353	to_be_resorted => ti
1354	ELSE
1355	IF ( .NOT. ALLOCATED( ti_av ) ) THEN
1356	ALLOCATE( ti_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1357	ti_av = REAL( fill_value, KIND = wp )
1358	ENDIF
1359	to_be_resorted => ti_av
1360	ENDIF
1361	flag_nr = 0
1362
1363	CASE ( 'uu' )
1364	IF ( av == 0 ) THEN
1365	to_be_resorted => uu
1366	ELSE
1367	IF ( .NOT. ALLOCATED( uu_av ) ) THEN
1368	ALLOCATE( uu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1369	uu_av = REAL( fill_value, KIND = wp )
1370	ENDIF
1371	to_be_resorted => uu_av
1372	ENDIF
1373	flag_nr = 1
1374
1375	CASE ( 'vv' )
1376	IF ( av == 0 ) THEN
1377	to_be_resorted => vv
1378	ELSE
1379	IF ( .NOT. ALLOCATED( vv_av ) ) THEN
1380	ALLOCATE( vv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1381	vv_av = REAL( fill_value, KIND = wp )
1382	ENDIF
1383	to_be_resorted => vv_av
1384	ENDIF
1385	flag_nr = 2
1386
1387	CASE ( 'ww' )
1388	IF ( av == 0 ) THEN
1389	to_be_resorted => ww
1390	ELSE
1391	IF ( .NOT. ALLOCATED( ww_av ) ) THEN
1392	ALLOCATE( ww_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1393	ww_av = REAL( fill_value, KIND = wp )
1394	ENDIF
1395	to_be_resorted => ww_av
1396	ENDIF
1397	flag_nr = 3
1398
1399	CASE ( 'wu' )
1400	IF ( av == 0 ) THEN
1401	to_be_resorted => wu
1402	ELSE
1403	IF ( .NOT. ALLOCATED( wu_av ) ) THEN
1404	ALLOCATE( wu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1405	wu_av = REAL( fill_value, KIND = wp )
1406	ENDIF
1407	to_be_resorted => wu_av
1408	ENDIF
1409	flag_nr = 0
1410
1411	CASE ( 'wv' )
1412	IF ( av == 0 ) THEN
1413	to_be_resorted => wv
1414	ELSE
1415	IF ( .NOT. ALLOCATED( wv_av ) ) THEN
1416	ALLOCATE( wv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1417	wv_av = REAL( fill_value, KIND = wp )
1418	ENDIF
1419	to_be_resorted => wv_av
1420	ENDIF
1421	flag_nr = 0
1422
1423	CASE ( 'wtheta' )
1424	IF ( av == 0 ) THEN
1425	to_be_resorted => wtheta
1426	ELSE
1427	IF ( .NOT. ALLOCATED( wtheta_av ) ) THEN
1428	ALLOCATE( wtheta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1429	wtheta_av = REAL( fill_value, KIND = wp )
1430	ENDIF
1431	to_be_resorted => wtheta_av
1432	ENDIF
1433	flag_nr = 0
1434
1435	CASE ( 'wq' )
1436	IF ( av == 0 ) THEN
1437	to_be_resorted => wq
1438	ELSE
1439	IF ( .NOT. ALLOCATED( wq_av ) ) THEN
1440	ALLOCATE( wq_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1441	wq_av = REAL( fill_value, KIND = wp )
1442	ENDIF
1443	to_be_resorted => wq_av
1444	ENDIF
1445	flag_nr = 0
1446
1447	CASE ( 'ws' )
1448	IF ( av == 0 ) THEN
1449	to_be_resorted => ws
1450	ELSE
1451	IF ( .NOT. ALLOCATED( ws_av ) ) THEN
1452	ALLOCATE( ws_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1453	ws_av = REAL( fill_value, KIND = wp )
1454	ENDIF
1455	to_be_resorted => ws_av
1456	ENDIF
1457	flag_nr = 0
1458
1459	CASE ( 'wspeed' )
1460	IF ( av == 0 ) THEN
1461	to_be_resorted => wspeed
1462	ELSE
1463	IF ( .NOT. ALLOCATED( wspeed_av ) ) THEN
1464	ALLOCATE( wspeed_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1465	wspeed_av = REAL( fill_value, KIND = wp )
1466	ENDIF
1467	to_be_resorted => wspeed_av
1468	ENDIF
1469	flag_nr = 0
1470
1471	CASE ( 'wdir' )
1472	IF ( av == 0 ) THEN
1473	to_be_resorted => wdir
1474	ELSE
1475	IF ( .NOT. ALLOCATED( wdir_av ) ) THEN
1476	ALLOCATE( wdir_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1477	wdir_av = REAL( fill_value, KIND = wp )
1478	ENDIF
1479	to_be_resorted => wdir_av
1480	ENDIF
1481	flag_nr = 0
1482
1483	CASE DEFAULT
1484	found = .FALSE.
1485
1486	END SELECT
1487
1488	IF ( found .AND. .NOT. resorted ) THEN
1489	DO i = nxl, nxr
1490	DO j = nys, nyn
1491	DO k = nzb_do, nzt_do
1492	local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), &
1493	REAL( fill_value, KIND = wp ), &
1494	BTEST( wall_flags_total_0(k,j,i), flag_nr ) )
1495	ENDDO
1496	ENDDO
1497	ENDDO
1498	ENDIF
1499
1500	END SUBROUTINE doq_output_3d
1501
1502
1503	!--------------------------------------------------------------------------------------------------!
1504	!
1505	! Description:
1506	! ------------
1507	!> Resorts the user-defined output quantity with indices (k,j,i) to a temporary array with indices
1508	!> (i,j,k) for masked data output.
1509	!--------------------------------------------------------------------------------------------------!
1510	SUBROUTINE doq_output_mask( av, variable, found, local_pf, mid )
1511
1512	USE control_parameters
1513
1514	USE indices
1515
1516	IMPLICIT NONE
1517
1518	REAL(wp), PARAMETER :: fill_value = -9999.0_wp !< value for the _FillValue attribute
1519
1520	CHARACTER (LEN=*) :: variable !<
1521	CHARACTER (LEN=5) :: grid !< flag to distinquish between staggered grids
1522
1523	INTEGER(iwp) :: av !< index indicating averaged or instantaneous output
1524	INTEGER(iwp) :: flag_nr !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
1525	INTEGER(iwp) :: i !< index variable along x-direction
1526	INTEGER(iwp) :: j !< index variable along y-direction
1527	INTEGER(iwp) :: k !< index variable along z-direction
1528	INTEGER(iwp) :: im !< loop index for masked variables
1529	INTEGER(iwp) :: jm !< loop index for masked variables
1530	INTEGER(iwp) :: kk !< masked output index variable along z-direction
1531	INTEGER(iwp) :: mid !< masked output running index
1532	INTEGER(iwp) :: ktt !< k index of highest horizontal surface
1533
1534	LOGICAL :: found !< true if variable is in list
1535	LOGICAL :: resorted !< true if array is resorted
1536
1537	REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) :: local_pf !<
1538
1539	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which needs to be resorted for output
1540
1541
1542	flag_nr = 0
1543	found = .TRUE.
1544	resorted = .FALSE.
1545	grid = 's'
1546
1547	SELECT CASE ( TRIM( variable ) )
1548
1549	CASE ( 'rh' )
1550	IF ( av == 0 ) THEN
1551	to_be_resorted => rh
1552	ELSE
1553	to_be_resorted => rh_av
1554	ENDIF
1555	grid = 's'
1556	flag_nr = 0
1557
1558	CASE ( 'ta' )
1559	IF ( av == 0 ) THEN
1560	to_be_resorted => ta
1561	ELSE
1562	to_be_resorted => ta_av
1563	ENDIF
1564	grid = 's'
1565	flag_nr = 0
1566
1567	CASE ( 'ti' )
1568	IF ( av == 0 ) THEN
1569	to_be_resorted => ti
1570	ELSE
1571	to_be_resorted => ti_av
1572	ENDIF
1573	grid = 's'
1574	flag_nr = 0
1575
1576	CASE ( 'uu' )
1577	IF ( av == 0 ) THEN
1578	to_be_resorted => uu
1579	ELSE
1580	to_be_resorted => uu_av
1581	ENDIF
1582	grid = 'u'
1583	flag_nr = 1
1584
1585	CASE ( 'vv' )
1586	IF ( av == 0 ) THEN
1587	to_be_resorted => vv
1588	ELSE
1589	to_be_resorted => vv_av
1590	ENDIF
1591	grid = 'v'
1592	flag_nr = 2
1593
1594	CASE ( 'ww' )
1595	IF ( av == 0 ) THEN
1596	to_be_resorted => ww
1597	ELSE
1598	to_be_resorted => ww_av
1599	ENDIF
1600	grid = 'w'
1601	flag_nr = 3
1602
1603	CASE ( 'wu' )
1604	IF ( av == 0 ) THEN
1605	to_be_resorted => wu
1606	ELSE
1607	to_be_resorted => wu_av
1608	ENDIF
1609	grid = 's'
1610	flag_nr = 0
1611
1612	CASE ( 'wv' )
1613	IF ( av == 0 ) THEN
1614	to_be_resorted => wv
1615	ELSE
1616	to_be_resorted => wv_av
1617	ENDIF
1618	grid = 's'
1619	flag_nr = 0
1620
1621	CASE ( 'wtheta' )
1622	IF ( av == 0 ) THEN
1623	to_be_resorted => wtheta
1624	ELSE
1625	to_be_resorted => wtheta_av
1626	ENDIF
1627	grid = 's'
1628	flag_nr = 0
1629
1630	CASE ( 'wq' )
1631	IF ( av == 0 ) THEN
1632	to_be_resorted => wq
1633	ELSE
1634	to_be_resorted => wq_av
1635	ENDIF
1636	grid = 's'
1637	flag_nr = 0
1638
1639	CASE ( 'ws' )
1640	IF ( av == 0 ) THEN
1641	to_be_resorted => ws
1642	ELSE
1643	to_be_resorted => ws_av
1644	ENDIF
1645	grid = 's'
1646	flag_nr = 0
1647
1648	CASE ( 'wspeed' )
1649	IF ( av == 0 ) THEN
1650	to_be_resorted => wspeed
1651	ELSE
1652	to_be_resorted => wspeed_av
1653	ENDIF
1654	grid = 's'
1655	flag_nr = 0
1656
1657	CASE ( 'wdir' )
1658	IF ( av == 0 ) THEN
1659	to_be_resorted => wdir
1660	ELSE
1661	to_be_resorted => wdir_av
1662	ENDIF
1663	grid = 's'
1664	flag_nr = 0
1665
1666	CASE DEFAULT
1667	found = .FALSE.
1668
1669	END SELECT
1670
1671	IF ( found .AND. .NOT. resorted ) THEN
1672	IF ( .NOT. mask_surface(mid) ) THEN
1673	!
1674	!-- Default masked output
1675	DO i = 1, mask_size_l(mid,1)
1676	DO j = 1, mask_size_l(mid,2)
1677	DO k = 1, mask_size_l(mid,3)
1678	local_pf(i,j,k) = MERGE( to_be_resorted(mask_k(mid,k), &
1679	mask_j(mid,j), &
1680	mask_i(mid,i)), &
1681	REAL( fill_value, KIND = wp ), &
1682	BTEST( wall_flags_total_0(mask_k(mid,k), &
1683	mask_j(mid,j), &
1684	mask_i(mid,i)), &
1685	flag_nr ) )
1686	ENDDO
1687	ENDDO
1688	ENDDO
1689
1690	ELSE
1691	!
1692	!-- Terrain-following masked output
1693	DO i = 1, mask_size_l(mid,1)
1694	DO j = 1, mask_size_l(mid,2)
1695	!
1696	!-- Get k index of the lowest non-terrain grid point, else set fill value
1697	im = mask_i(mid,i)
1698	jm = mask_j(mid,j)
1699	ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 ) ), DIM=1 ) - 1
1700	DO k = 1, mask_size_l(mid,3)
1701	kk = MIN( ktt + mask_k(mid,k) - 1, nzt+1 )
1702	!
1703	!-- Set value if not in building
1704	IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) ) THEN
1705	local_pf(i,j,k) = fill_value
1706	ELSE
1707	local_pf(i,j,k) = to_be_resorted(kk,jm,im)
1708	ENDIF
1709	ENDDO
1710	ENDDO
1711	ENDDO
1712
1713	ENDIF
1714	ENDIF
1715
1716	END SUBROUTINE doq_output_mask
1717
1718	!--------------------------------------------------------------------------------------------------!
1719	! Description:
1720	! ------------
1721	!> Allocate required arrays
1722	!--------------------------------------------------------------------------------------------------!
1723	SUBROUTINE doq_init
1724
1725	IMPLICIT NONE
1726
1727	INTEGER(iwp) :: ivar !< loop index over all 2d/3d/mask output quantities
1728
1729
1730	!
1731	!-- Next line is to avoid compiler warnings about unused variables
1732	IF ( timestep_number_at_prev_calc == 0 ) CONTINUE
1733	!
1734	!-- Preparatory steps and initialization of output arrays
1735	IF ( .NOT. prepared_diagnostic_output_quantities ) CALL doq_prepare
1736
1737	initialized_diagnostic_output_quantities = .FALSE.
1738
1739	ivar = 1
1740
1741	DO WHILE ( ivar <= SIZE( do_all ) )
1742
1743	SELECT CASE ( TRIM( do_all(ivar) ) )
1744	!
1745	!-- Allocate array for 'relative humidity'
1746	CASE ( 'rh' )
1747	IF ( .NOT. ALLOCATED( rh ) ) THEN
1748	ALLOCATE( rh(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1749	rh = 0.0_wp
1750	ENDIF
1751	!
1752	!-- Allocate array for 'air temperature'
1753	CASE ( 'ta' )
1754	IF ( .NOT. ALLOCATED( ta ) ) THEN
1755	ALLOCATE( ta(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1756	ta = 0.0_wp
1757	ENDIF
1758	!
1759	!-- Allocate array for 'turbulence intensity'
1760	CASE ( 'ti' )
1761	IF ( .NOT. ALLOCATED( ti ) ) THEN
1762	ALLOCATE( ti(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1763	ti = 0.0_wp
1764	ENDIF
1765	!
1766	!-- Allocate array for uu
1767	CASE ( 'uu' )
1768	IF ( .NOT. ALLOCATED( uu ) ) THEN
1769	ALLOCATE( uu(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1770	uu = 0.0_wp
1771	ENDIF
1772	!
1773	!-- Allocate array for vv
1774	CASE ( 'vv' )
1775	IF ( .NOT. ALLOCATED( vv ) ) THEN
1776	ALLOCATE( vv(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1777	vv = 0.0_wp
1778	ENDIF
1779	!
1780	!-- Allocate array for ww
1781	CASE ( 'ww' )
1782	IF ( .NOT. ALLOCATED( ww ) ) THEN
1783	ALLOCATE( ww(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1784	ww = 0.0_wp
1785	ENDIF
1786	!
1787	!-- Allocate array for wu
1788	CASE ( 'wu' )
1789	IF ( .NOT. ALLOCATED( wu ) ) THEN
1790	ALLOCATE( wu(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1791	wu = 0.0_wp
1792	ENDIF
1793	!
1794	!-- Allocate array for wv
1795	CASE ( 'wv' )
1796	IF ( .NOT. ALLOCATED( wv ) ) THEN
1797	ALLOCATE( wv(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1798	wv = 0.0_wp
1799	ENDIF
1800	!
1801	!-- Allocate array for wtheta
1802	CASE ( 'wtheta' )
1803	IF ( .NOT. ALLOCATED( wtheta ) ) THEN
1804	ALLOCATE( wtheta(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1805	wtheta = 0.0_wp
1806	ENDIF
1807	!
1808	!-- Allocate array for wq
1809	CASE ( 'wq' )
1810	IF ( .NOT. ALLOCATED( wq ) ) THEN
1811	ALLOCATE( wq(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1812	wq = 0.0_wp
1813	ENDIF
1814	!
1815	!-- Allocate array for ws
1816	CASE ( 'ws' )
1817	IF ( .NOT. ALLOCATED( ws ) ) THEN
1818	ALLOCATE( ws(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1819	ws = 0.0_wp
1820	ENDIF
1821	!
1822	!-- Allocate array for 2-m potential temperature
1823	CASE ( 'theta_2m*' )
1824	IF ( .NOT. ALLOCATED( pt_2m ) ) THEN
1825	ALLOCATE( pt_2m(nysg:nyng,nxlg:nxrg) )
1826	pt_2m = 0.0_wp
1827	ENDIF
1828	!
1829	!-- Allocate array for 10-m wind speed
1830	CASE ( 'wspeed_10m*' )
1831	IF ( .NOT. ALLOCATED( uv_10m ) ) THEN
1832	ALLOCATE( uv_10m(nysg:nyng,nxlg:nxrg) )
1833	uv_10m = 0.0_wp
1834	ENDIF
1835	!
1836	!-- Allocate array for wspeed
1837	CASE ( 'wspeed' )
1838	IF ( .NOT. ALLOCATED( wspeed ) ) THEN
1839	ALLOCATE( wspeed(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1840	wspeed = 0.0_wp
1841	ENDIF
1842
1843	!
1844	!-- Allocate array for wdir
1845	CASE ( 'wdir' )
1846	IF ( .NOT. ALLOCATED( u_center ) ) THEN
1847	ALLOCATE( u_center(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1848	u_center = 0.0_wp
1849	ENDIF
1850	IF ( .NOT. ALLOCATED( v_center ) ) THEN
1851	ALLOCATE( v_center(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1852	v_center = 0.0_wp
1853	ENDIF
1854	IF ( .NOT. ALLOCATED( wdir ) ) THEN
1855	ALLOCATE( wdir(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1856	wdir = 0.0_wp
1857	ENDIF
1858
1859	END SELECT
1860
1861	ivar = ivar + 1
1862	ENDDO
1863
1864	initialized_diagnostic_output_quantities = .TRUE.
1865
1866	END SUBROUTINE doq_init
1867
1868
1869	!--------------------------------------------------------------------------------------------------!
1870	! Description:
1871	! ------------
1872	!> Calculation of module-specific statistics, i.e. horizontally averaged profiles and time series.
1873	!> This is called for every statistic region sr, but at least for the region "total domain" (sr=0).
1874	!--------------------------------------------------------------------------------------------------!
1875	SUBROUTINE doq_statistics( mode, sr, tn )
1876
1877
1878	CHARACTER (LEN=*) :: mode !<
1879
1880	! INTEGER(iwp) :: i !<
1881	! INTEGER(iwp) :: j !<
1882	! INTEGER(iwp) :: k !<
1883	INTEGER(iwp) :: sr !<
1884	INTEGER(iwp) :: tn !<
1885	!
1886	!-- Next line is to avoid compiler warning about unused variables. Please remove.
1887	IF ( sr == 0 .OR. tn == 0 ) CONTINUE
1888
1889	IF ( mode == 'profiles' ) THEN
1890
1891	ELSEIF ( mode == 'time_series' ) THEN
1892
1893	ENDIF
1894
1895	END SUBROUTINE doq_statistics
1896
1897
1898	!--------------------------------------------------------------------------------------------------!
1899	! Description:
1900	! ------------
1901	!> Calculate diagnostic quantities
1902	!--------------------------------------------------------------------------------------------------!
1903	SUBROUTINE doq_calculate
1904
1905	IMPLICIT NONE
1906
1907	INTEGER(iwp) :: i !< grid index x-dimension
1908	INTEGER(iwp) :: j !< grid index y-dimension
1909	INTEGER(iwp) :: k !< grid index z-dimension
1910	INTEGER(iwp) :: ivar !< loop index over all 2d/3d/mask output quantities
1911
1912	REAL(wp) :: temp !< temperature
1913	REAL(wp) :: e_s !< saturation vapor pressure
1914	REAL(wp) :: q_s !< saturation mixing ratio
1915
1916	TYPE(surf_type), POINTER :: surf !< surf-type array, used to generalize subroutines
1917
1918
1919	! CALL cpu_log( log_point(41), 'calculate_quantities', 'start' )
1920
1921	!
1922	!-- Save timestep number to check in time_integration if doq_calculate has been called already,
1923	!-- since the CALL occurs at two locations, but the calculations need to be done only once per
1924	!-- timestep.
1925	timestep_number_at_prev_calc = current_timestep_number
1926
1927	ivar = 1
1928
1929	DO WHILE ( ivar <= SIZE( do_all ) )
1930
1931	SELECT CASE ( TRIM( do_all(ivar) ) )
1932	!
1933	!-- rh (relative humidity)
1934	CASE ( 'rh' )
1935	DO i = nxl, nxr
1936	DO j = nys, nyn
1937	DO k = nzb+1, nzt
1938	IF ( humidity ) THEN
1939	temp = exner(k) * pt(k,j,i)
1940	e_s = magnus( temp )
1941	q_s = rd_d_rv * e_s / ( hyp(k) - e_s )
1942	rh(k,j,i) = q(k,j,i) / q_s * 100.0_wp &
1943	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
1944	ENDIF
1945	ENDDO
1946	ENDDO
1947	ENDDO
1948	!
1949	!-- ta (air temperature)
1950	CASE ( 'ta' )
1951	DO i = nxl, nxr
1952	DO j = nys, nyn
1953	DO k = nzb+1, nzt
1954	ta(k,j,i) = exner(k) * pt(k,j,i) - degc_to_k
1955	ENDDO
1956	ENDDO
1957	ENDDO
1958	!
1959	!-- Calculate 'turbulence intensity' from rot[(u,v,w)] at scalar grid point
1960	CASE ( 'ti' )
1961	DO i = nxl, nxr
1962	DO j = nys, nyn
1963	DO k = nzb+1, nzt
1964	ti(k,j,i) = 0.25_wp * SQRT( &
1965	( ( w(k,j+1,i) + w(k-1,j+1,i) &
1966	- w(k,j-1,i) - w(k-1,j-1,i) ) * ddy &
1967	- ( v(k+1,j,i) + v(k+1,j+1,i) &
1968	- v(k-1,j,i) - v(k-1,j+1,i) ) * ddzu(k) )**2 &
1969	+ ( ( u(k+1,j,i) + u(k+1,j,i+1) &
1970	- u(k-1,j,i) - u(k-1,j,i+1) ) * ddzu(k) &
1971	- ( w(k,j,i+1) + w(k-1,j,i+1) &
1972	- w(k,j,i-1) - w(k-1,j,i-1) ) * ddx )**2 &
1973	+ ( ( v(k,j,i+1) + v(k,j+1,i+1) &
1974	- v(k,j,i-1) - v(k,j+1,i-1) ) * ddx &
1975	- ( u(k,j+1,i) + u(k,j+1,i+1) &
1976	- u(k,j-1,i) - u(k,j-1,i+1) ) * ddy )**2 ) &
1977	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
1978	ENDDO
1979	ENDDO
1980	ENDDO
1981	!
1982	!-- uu
1983	CASE ( 'uu' )
1984	DO i = nxl, nxr
1985	DO j = nys, nyn
1986	DO k = nzb+1, nzt
1987	uu(k,j,i) = u(k,j,i) * u(k,j,i) &
1988	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
1989	ENDDO
1990	ENDDO
1991	ENDDO
1992	!
1993	!-- vv
1994	CASE ( 'vv' )
1995	DO i = nxl, nxr
1996	DO j = nys, nyn
1997	DO k = nzb+1, nzt
1998	vv(k,j,i) = v(k,j,i) * v(k,j,i) &
1999	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
2000	ENDDO
2001	ENDDO
2002	ENDDO
2003	!
2004	!-- ww
2005	CASE ( 'ww' )
2006	DO i = nxl, nxr
2007	DO j = nys, nyn
2008	DO k = nzb+1, nzt-1
2009	ww(k,j,i) = w(k,j,i) * w(k,j,i) &
2010	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
2011	ENDDO
2012	ENDDO
2013	ENDDO
2014	!
2015	!-- wu
2016	CASE ( 'wu' )
2017	DO i = nxl, nxr
2018	DO j = nys, nyn
2019	DO k = nzb+1, nzt-1
2020	wu(k,j,i) = w(k,j,i) &
2021	* 0.25_wp * ( u(k,j,i) + u(k,j,i+1) + u(k+1,j,i) + u(k+1,j,i+1) )&
2022	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
2023	ENDDO
2024	ENDDO
2025	ENDDO
2026	!
2027	!-- wv
2028	CASE ( 'wv' )
2029	DO i = nxl, nxr
2030	DO j = nys, nyn
2031	DO k = nzb+1, nzt-1
2032	wv(k,j,i) = w(k,j,i) &
2033	* 0.25_wp * ( v(k,j,i) + v(k,j+1,i) + v(k+1,j,i) + v(k+1,j+1,i) )&
2034	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
2035	ENDDO
2036	ENDDO
2037	ENDDO
2038	!
2039	!-- wtheta
2040	CASE ( 'wtheta' )
2041	DO i = nxl, nxr
2042	DO j = nys, nyn
2043	DO k = nzb+1, nzt-1
2044	wtheta(k,j,i) = w(k,j,i) &
2045	* 0.5_wp * ( pt(k,j,i) + pt(k+1,j,i) ) &
2046	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ))
2047	ENDDO
2048	ENDDO
2049	ENDDO
2050	!
2051	!-- wq
2052	CASE ( 'wq' )
2053	DO i = nxl, nxr
2054	DO j = nys, nyn
2055	DO k = nzb+1, nzt-1
2056	wq(k,j,i) = w(k,j,i) * 0.5_wp * ( q(k,j,i) + q(k+1,j,i) ) &
2057	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
2058	ENDDO
2059	ENDDO
2060	ENDDO
2061	!
2062	!-- ws
2063	CASE ( 'ws' )
2064	DO i = nxl, nxr
2065	DO j = nys, nyn
2066	DO k = nzb+1, nzt-1
2067	ws(k,j,i) = w(k,j,i) * 0.5_wp * ( s(k,j,i) + s(k+1,j,i) ) &
2068	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
2069	ENDDO
2070	ENDDO
2071	ENDDO
2072	!
2073	!-- 2-m potential temperature
2074	CASE ( 'theta_2m*' )
2075	!
2076	!-- 2-m potential temperature is caluclated from surface arrays. In case the 2m level is
2077	!-- below the first grid point, MOST is applied, else, linear interpolation between two
2078	!-- vertical grid levels is applied. To access all surfaces, iterate over all horizontally-
2079	!-- upward facing surface types.
2080	surf => surf_def_h(0)
2081	CALL calc_pt_2m
2082	surf => surf_lsm_h(0)
2083	CALL calc_pt_2m
2084	surf => surf_usm_h(0)
2085	CALL calc_pt_2m
2086	!
2087	!-- 10-m wind speed
2088	CASE ( 'wspeed_10m*' )
2089	!
2090	!-- 10-m wind speed is caluclated from surface arrays. In case the 10m level is below the
2091	!-- first grid point, MOST is applied, else, linear interpolation between two vertical grid
2092	!-- levels is applied. To access all surfaces, iterate over all horizontally-upward facing
2093	!-- surface types.
2094	surf => surf_def_h(0)
2095	CALL calc_wind_10m
2096	surf => surf_lsm_h(0)
2097	CALL calc_wind_10m
2098	surf => surf_usm_h(0)
2099	CALL calc_wind_10m
2100	!
2101	!-- Horizontal wind speed
2102	CASE ( 'wspeed' )
2103	DO i = nxl, nxr
2104	DO j = nys, nyn
2105	DO k = nzb, nzt+1
2106	wspeed(k,j,i) = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 &
2107	+ ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2 ) &
2108	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ))
2109	ENDDO
2110	ENDDO
2111	ENDDO
2112
2113	!
2114	!-- Horizontal wind direction
2115	CASE ( 'wdir' )
2116	DO i = nxl, nxr
2117	DO j = nys, nyn
2118	DO k = nzb, nzt+1
2119	u_center(k,j,i) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
2120	v_center(k,j,i) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
2121
2122	wdir(k,j,i) = ATAN2( u_center(k,j,i), v_center(k,j,i) ) &
2123	/ pi * 180.0_wp + 180.0_wp
2124	ENDDO
2125	ENDDO
2126	ENDDO
2127
2128	END SELECT
2129
2130	ivar = ivar + 1
2131	ENDDO
2132
2133	! CALL cpu_log( log_point(41), 'calculate_quantities', 'stop' )
2134
2135	!
2136	!-- The following block contains subroutines to calculate diagnostic surface quantities.
2137	CONTAINS
2138	!--------------------------------------------------------------------------------------------------!
2139	! Description:
2140	! ------------
2141	!> Calculation of 2-m potential temperature.
2142	!--------------------------------------------------------------------------------------------------!
2143	SUBROUTINE calc_pt_2m
2144
2145	USE surface_layer_fluxes_mod, &
2146	ONLY: psi_h
2147
2148	IMPLICIT NONE
2149
2150	INTEGER(iwp) :: kk !< running index along the z-dimension
2151	INTEGER(iwp) :: m !< running index for surface elements
2152
2153	DO m = 1, surf%ns
2154
2155	i = surf%i(m)
2156	j = surf%j(m)
2157	k = surf%k(m)
2158	!
2159	!-- If 2-m level is below the first grid level, MOST is used for calculation of
2160	!-- 2-m temperature.
2161	IF ( surf%z_mo(m) > 2.0_wp ) THEN
2162	pt_2m(j,i) = surf%pt_surface(m) + surf%ts(m) / kappa &
2163	* ( LOG( 2.0_wp / surf%z0h(m) ) &
2164	- psi_h( 2.0_wp / surf%ol(m) ) &
2165	+ psi_h( surf%z0h(m) / surf%ol(m) ) )
2166	!
2167	!-- If 2-m level is above the first grid level, 2-m temperature is linearly interpolated
2168	!-- between the two nearest vertical grid levels. Note, since 2-m temperature is only
2169	!-- computed for horizontal upward-facing surfaces, only a vertical interpolation is
2170	!-- necessary.
2171	ELSE
2172	!
2173	!-- zw(k-1) defines the height of the surface.
2174	kk = k
2175	DO WHILE ( zu(kk) - zw(k-1) < 2.0_wp .AND. kk <= nzt )
2176	kk = kk + 1
2177	ENDDO
2178	!
2179	!-- kk defines the index of the first grid level >= 2m.
2180	pt_2m(j,i) = pt(kk-1,j,i) + &
2181	( zw(k-1) + 2.0_wp - zu(kk-1) ) * &
2182	( pt(kk,j,i) - pt(kk-1,j,i) ) / &
2183	( zu(kk) - zu(kk-1) )
2184	ENDIF
2185
2186	ENDDO
2187
2188	END SUBROUTINE calc_pt_2m
2189
2190	!--------------------------------------------------------------------------------------------------!
2191	! Description:
2192	! ------------
2193	!> Calculation of 10-m wind speed.
2194	!--------------------------------------------------------------------------------------------------!
2195	SUBROUTINE calc_wind_10m
2196
2197	USE surface_layer_fluxes_mod, &
2198	ONLY: psi_m
2199
2200	IMPLICIT NONE
2201
2202	INTEGER(iwp) :: kk !< running index along the z-dimension
2203	INTEGER(iwp) :: m !< running index for surface elements
2204
2205	REAL(wp) :: uv_l !< wind speed at lower grid point
2206	REAL(wp) :: uv_u !< wind speed at upper grid point
2207
2208
2209	DO m = 1, surf%ns
2210
2211	i = surf%i(m)
2212	j = surf%j(m)
2213	k = surf%k(m)
2214	!
2215	!-- If 10-m level is below the first grid level, MOST is used for calculation of 10-m
2216	!-- temperature.
2217	IF ( surf%z_mo(m) > 10.0_wp ) THEN
2218	uv_10m(j,i) = surf%us(m) / kappa &
2219	* ( LOG( 10.0_wp / surf%z0(m) ) &
2220	- psi_m( 10.0_wp / surf%ol(m) ) &
2221	+ psi_m( surf%z0(m) / surf%ol(m) ) )
2222	!
2223	!-- If 10-m level is above the first grid level, 10-m wind speed is linearly interpolated
2224	!-- between the two nearest vertical grid levels. Note, since 10-m temperature is only
2225	!-- computed for horizontal upward-facing surfaces, only a vertical interpolation is
2226	!-- necessary.
2227	ELSE
2228	!
2229	!-- zw(k-1) defines the height of the surface.
2230	kk = k
2231	DO WHILE ( zu(kk) - zw(k-1) < 10.0_wp .AND. kk <= nzt )
2232	kk = kk + 1
2233	ENDDO
2234	!
2235	!-- kk defines the index of the first grid level >= 10m.
2236	uv_l = SQRT( ( 0.5_wp * ( u(kk-1,j,i) + u(kk-1,j,i+1) ) )**2 &
2237	+ ( 0.5_wp * ( v(kk-1,j,i) + v(kk-1,j+1,i) ) )**2 )
2238
2239	uv_u = SQRT( ( 0.5_wp * ( u(kk,j,i) + u(kk,j,i+1) ) )**2 &
2240	+ ( 0.5_wp * ( v(kk,j,i) + v(kk,j+1,i) ) )**2 )
2241
2242	uv_10m(j,i) = uv_l + ( zw(k-1) + 10.0_wp - zu(kk-1) ) * &
2243	( uv_u - uv_l ) / &
2244	( zu(kk) - zu(kk-1) )
2245
2246	ENDIF
2247
2248	ENDDO
2249
2250	END SUBROUTINE calc_wind_10m
2251
2252	END SUBROUTINE doq_calculate
2253
2254
2255	!--------------------------------------------------------------------------------------------------!
2256	! Description:
2257	! ------------
2258	!> Preparation of the diagnostic output, counting of the module-specific output quantities and
2259	!> gathering of the output names.
2260	!--------------------------------------------------------------------------------------------------!
2261	SUBROUTINE doq_prepare
2262
2263	USE control_parameters, &
2264	ONLY: do2d, do3d, domask, masks
2265
2266	IMPLICIT NONE
2267
2268	CHARACTER (LEN=varnamelength), DIMENSION(0:1,500) :: do2d_var = ' ' !< label array for 2d output quantities
2269
2270	INTEGER(iwp) :: av !< index defining type of output, av=0 instantaneous, av=1 averaged
2271	INTEGER(iwp) :: ivar !< loop index
2272	INTEGER(iwp) :: ivar_all !< loop index
2273	INTEGER(iwp) :: l !< index for cutting string
2274	INTEGER(iwp) :: mid !< masked output running index
2275
2276
2277	prepared_diagnostic_output_quantities = .FALSE.
2278
2279	ivar = 1
2280	ivar_all = 1
2281
2282	DO av = 0, 1
2283	!
2284	!-- Remove _xy, _xz, or _yz from string
2285	l = MAX( 3, LEN_TRIM( do2d(av,ivar) ) )
2286	do2d_var(av,ivar)(1:l-3) = do2d(av,ivar)(1:l-3)
2287	!
2288	!-- Gather 2d output quantity names.
2289	!-- Check for double occurrence of output quantity, e.g. by _xy, _yz, _xz.
2290	DO WHILE ( do2d_var(av,ivar)(1:1) /= ' ' )
2291	IF ( .NOT. ANY( do_all == do2d_var(av,ivar) ) ) THEN
2292	do_all(ivar_all) = do2d_var(av,ivar)
2293	ENDIF
2294	ivar = ivar + 1
2295	ivar_all = ivar_all + 1
2296	l = MAX( 3, LEN_TRIM( do2d(av,ivar) ) )
2297	do2d_var(av,ivar)(1:l-3) = do2d(av,ivar)(1:l-3)
2298	ENDDO
2299
2300	ivar = 1
2301	!
2302	!-- Gather 3d output quantity names
2303	DO WHILE ( do3d(av,ivar)(1:1) /= ' ' )
2304	do_all(ivar_all) = do3d(av,ivar)
2305	ivar = ivar + 1
2306	ivar_all = ivar_all + 1
2307	ENDDO
2308
2309	ivar = 1
2310	!
2311	!-- Gather masked output quantity names. Also check for double output e.g. by different masks.
2312	DO mid = 1, masks
2313	DO WHILE ( domask(mid,av,ivar)(1:1) /= ' ' )
2314	IF ( .NOT. ANY( do_all == domask(mid,av,ivar) ) ) THEN
2315	do_all(ivar_all) = domask(mid,av,ivar)
2316	ENDIF
2317
2318	ivar = ivar + 1
2319	ivar_all = ivar_all + 1
2320	ENDDO
2321	ivar = 1
2322	ENDDO
2323
2324	ENDDO
2325
2326	prepared_diagnostic_output_quantities = .TRUE.
2327
2328	END SUBROUTINE doq_prepare
2329
2330	!--------------------------------------------------------------------------------------------------!
2331	! Description:
2332	! ------------
2333	!> Subroutine reads local (subdomain) restart data
2334	!--------------------------------------------------------------------------------------------------!
2335	SUBROUTINE doq_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, nxr_on_file, nynf, nync, &
2336	nyn_on_file, nysf, nysc, nys_on_file, tmp_2d, tmp_3d, found )
2337
2338	USE control_parameters
2339
2340	USE indices
2341
2342	USE kinds
2343
2344	USE pegrid
2345
2346
2347	IMPLICIT NONE
2348
2349	INTEGER(iwp) :: k !<
2350	INTEGER(iwp) :: nxlc !<
2351	INTEGER(iwp) :: nxlf !<
2352	INTEGER(iwp) :: nxl_on_file !<
2353	INTEGER(iwp) :: nxrc !<
2354	INTEGER(iwp) :: nxrf !<
2355	INTEGER(iwp) :: nxr_on_file !<
2356	INTEGER(iwp) :: nync !<
2357	INTEGER(iwp) :: nynf !<
2358	INTEGER(iwp) :: nyn_on_file !<
2359	INTEGER(iwp) :: nysc !<
2360	INTEGER(iwp) :: nysf !<
2361	INTEGER(iwp) :: nys_on_file !<
2362
2363	LOGICAL, INTENT(OUT) :: found
2364
2365	REAL(wp), DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_2d !<
2366
2367	REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d !<
2368
2369
2370	found = .TRUE.
2371
2372	SELECT CASE ( restart_string(1:length) )
2373
2374	CASE ( 'pt_2m_av' )
2375	IF ( .NOT. ALLOCATED( pt_2m_av ) ) THEN
2376	ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )
2377	ENDIF
2378	IF ( k == 1 ) READ ( 13 ) tmp_2d
2379	pt_2m_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
2380	tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2381
2382	CASE ( 'rh_av' )
2383	IF ( .NOT. ALLOCATED( rh_av ) ) THEN
2384	ALLOCATE( rh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2385	ENDIF
2386	IF ( k == 1 ) READ ( 13 ) tmp_3d
2387	rh_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
2388	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2389
2390	CASE ( 'ta_av' )
2391	IF ( .NOT. ALLOCATED( ta_av ) ) THEN
2392	ALLOCATE( ta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2393	ENDIF
2394	IF ( k == 1 ) READ ( 13 ) tmp_3d
2395	ta_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
2396	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2397
2398	CASE ( 'ti_av' )
2399	IF ( .NOT. ALLOCATED( ti_av ) ) THEN
2400	ALLOCATE( ti_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2401	ENDIF
2402	IF ( k == 1 ) READ ( 13 ) tmp_3d
2403	ti_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
2404	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2405
2406	CASE ( 'u_center_av' )
2407	IF ( .NOT. ALLOCATED( u_center_av ) ) THEN
2408	ALLOCATE( u_center_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2409	ENDIF
2410	IF ( k == 1 ) READ ( 13 ) tmp_3d
2411	u_center_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
2412	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2413
2414	CASE ( 'uu_av' )
2415	IF ( .NOT. ALLOCATED( uu_av ) ) THEN
2416	ALLOCATE( uu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2417	ENDIF
2418	IF ( k == 1 ) READ ( 13 ) tmp_3d
2419	uu_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
2420	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2421
2422	CASE ( 'uv_10m_av' )
2423	IF ( .NOT. ALLOCATED( uv_10m_av ) ) THEN
2424	ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) )
2425	ENDIF
2426	IF ( k == 1 ) READ ( 13 ) tmp_2d
2427	uv_10m_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
2428	tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2429
2430	CASE ( 'v_center_av' )
2431	IF ( .NOT. ALLOCATED( v_center_av ) ) THEN
2432	ALLOCATE( v_center_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2433	ENDIF
2434	IF ( k == 1 ) READ ( 13 ) tmp_3d
2435	v_center_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
2436	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2437
2438	CASE ( 'vv_av' )
2439	IF ( .NOT. ALLOCATED( vv_av ) ) THEN
2440	ALLOCATE( vv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2441	ENDIF
2442	IF ( k == 1 ) READ ( 13 ) tmp_3d
2443	vv_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
2444	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2445
2446	CASE ( 'wtheta_av' )
2447	IF ( .NOT. ALLOCATED( wtheta_av ) ) THEN
2448	ALLOCATE( wtheta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2449	ENDIF
2450	IF ( k == 1 ) READ ( 13 ) tmp_3d
2451	wtheta_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
2452	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2453
2454	CASE ( 'wq_av' )
2455	IF ( .NOT. ALLOCATED( wq_av ) ) THEN
2456	ALLOCATE( wq_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2457	ENDIF
2458	IF ( k == 1 ) READ ( 13 ) tmp_3d
2459	wq_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
2460	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2461
2462	CASE ( 'ws_av' )
2463	IF ( .NOT. ALLOCATED( ws_av ) ) THEN
2464	ALLOCATE( ws_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2465	ENDIF
2466	IF ( k == 1 ) READ ( 13 ) tmp_3d
2467	ws_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
2468	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2469
2470	CASE ( 'wspeed_av' )
2471	IF ( .NOT. ALLOCATED( wspeed_av ) ) THEN
2472	ALLOCATE( wspeed_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2473	ENDIF
2474	IF ( k == 1 ) READ ( 13 ) tmp_3d
2475	wspeed_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
2476	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2477
2478	CASE ( 'wu_av' )
2479	IF ( .NOT. ALLOCATED( wu_av ) ) THEN
2480	ALLOCATE( wu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2481	ENDIF
2482	IF ( k == 1 ) READ ( 13 ) tmp_3d
2483	wu_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
2484	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2485
2486	CASE ( 'wv_av' )
2487	IF ( .NOT. ALLOCATED( wv_av ) ) THEN
2488	ALLOCATE( wv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2489	ENDIF
2490	IF ( k == 1 ) READ ( 13 ) tmp_3d
2491	wv_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
2492	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2493
2494	CASE ( 'ww_av' )
2495	IF ( .NOT. ALLOCATED( ww_av ) ) THEN
2496	ALLOCATE( ww_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2497	ENDIF
2498	IF ( k == 1 ) READ ( 13 ) tmp_3d
2499	ww_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
2500	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
2501
2502
2503	CASE DEFAULT
2504
2505	found = .FALSE.
2506
2507	END SELECT
2508
2509	END SUBROUTINE doq_rrd_local_ftn
2510
2511	!--------------------------------------------------------------------------------------------------!
2512	! Description:
2513	! ------------
2514	!> Read module-specific local restart data arrays (MPI-IO).
2515	!--------------------------------------------------------------------------------------------------!
2516	SUBROUTINE doq_rrd_local_mpi
2517
2518	LOGICAL :: array_found !< control flad indicating whether the array is found in the file or not
2519
2520
2521	CALL rd_mpi_io_check_array( 'pt_2m_av' , found = array_found )
2522	IF ( array_found ) THEN
2523	IF ( .NOT. ALLOCATED( pt_2m_av ) ) ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )
2524	CALL rrd_mpi_io( 'pt_2m_av', pt_2m_av )
2525	ENDIF
2526
2527	CALL rd_mpi_io_check_array( 'ta_av' , found = array_found )
2528	IF ( array_found ) THEN
2529	IF ( .NOT. ALLOCATED( ta_av ) ) ALLOCATE( ta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2530	CALL rrd_mpi_io( 'ta_av', ta_av )
2531	ENDIF
2532
2533	CALL rd_mpi_io_check_array( 'ti_av' , found = array_found )
2534	IF ( array_found ) THEN
2535	IF ( .NOT. ALLOCATED( ti_av ) ) ALLOCATE( ti_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2536	CALL rrd_mpi_io( 'ti_av', ti_av )
2537	ENDIF
2538
2539	CALL rd_mpi_io_check_array( 'rh_av' , found = array_found )
2540	IF ( array_found ) THEN
2541	IF ( .NOT. ALLOCATED( rh_av ) ) ALLOCATE( ti_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2542	CALL rrd_mpi_io( 'rh_av', rh_av )
2543	ENDIF
2544
2545	CALL rd_mpi_io_check_array( 'u_center_av' , found = array_found )
2546	IF ( array_found ) THEN
2547	IF ( .NOT. ALLOCATED( u_center_av ) ) ALLOCATE( u_center_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2548	CALL rrd_mpi_io( 'u_center_av', u_center_av )
2549	ENDIF
2550
2551	CALL rd_mpi_io_check_array( 'uu_av' , found = array_found )
2552	IF ( array_found ) THEN
2553	IF ( .NOT. ALLOCATED( uu_av ) ) ALLOCATE( uu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2554	CALL rrd_mpi_io( 'uu_av', uu_av )
2555	ENDIF
2556
2557	CALL rd_mpi_io_check_array( 'uv_10m_av' , found = array_found )
2558	IF ( array_found ) THEN
2559	IF ( .NOT. ALLOCATED( uv_10m_av ) ) ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) )
2560	CALL rrd_mpi_io( 'uv_10m_av', uv_10m_av )
2561	ENDIF
2562
2563	CALL rd_mpi_io_check_array( 'v_center_av' , found = array_found )
2564	IF ( array_found ) THEN
2565	IF ( .NOT. ALLOCATED( v_center_av ) ) ALLOCATE( v_center_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2566	CALL rrd_mpi_io( 'v_center_av', v_center_av )
2567	ENDIF
2568
2569	CALL rd_mpi_io_check_array( 'vv_av' , found = array_found )
2570	IF ( array_found ) THEN
2571	IF ( .NOT. ALLOCATED( vv_av ) ) ALLOCATE( vv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2572	CALL rrd_mpi_io( 'vv_av', vv_av )
2573	ENDIF
2574
2575	CALL rd_mpi_io_check_array( 'ww_av' , found = array_found )
2576	IF ( array_found ) THEN
2577	IF ( .NOT. ALLOCATED( ww_av ) ) ALLOCATE( ww_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2578	CALL rrd_mpi_io( 'ww_av', ww_av )
2579	ENDIF
2580
2581	CALL rd_mpi_io_check_array( 'wu_av' , found = array_found )
2582	IF ( array_found ) THEN
2583	IF ( .NOT. ALLOCATED( wu_av ) ) ALLOCATE( wu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2584	CALL rrd_mpi_io( 'wu_av', wu_av )
2585	ENDIF
2586
2587	CALL rd_mpi_io_check_array( 'wv_av' , found = array_found )
2588	IF ( array_found ) THEN
2589	IF ( .NOT. ALLOCATED( wv_av ) ) ALLOCATE( wv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2590	CALL rrd_mpi_io( 'wv_av', wv_av )
2591	ENDIF
2592
2593	CALL rd_mpi_io_check_array( 'wtheta_av' , found = array_found )
2594	IF ( array_found ) THEN
2595	IF ( .NOT. ALLOCATED( wtheta_av ) ) ALLOCATE( wtheta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2596	CALL rrd_mpi_io( 'wtheta_av', wtheta_av )
2597	ENDIF
2598
2599	CALL rd_mpi_io_check_array( 'wq_av' , found = array_found )
2600	IF ( array_found ) THEN
2601	IF ( .NOT. ALLOCATED( wq_av ) ) ALLOCATE( wq_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2602	CALL rrd_mpi_io( 'wq_av', wq_av )
2603	ENDIF
2604
2605	CALL rd_mpi_io_check_array( 'ws_av' , found = array_found )
2606	IF ( array_found ) THEN
2607	IF ( .NOT. ALLOCATED( ws_av ) ) ALLOCATE( ws_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2608	CALL rrd_mpi_io( 'ws_av', ws_av )
2609	ENDIF
2610
2611	CALL rd_mpi_io_check_array( 'wspeed_av' , found = array_found )
2612	IF ( array_found ) THEN
2613	IF ( .NOT. ALLOCATED( wspeed_av ) ) ALLOCATE( wspeed_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
2614	CALL rrd_mpi_io( 'wspeed_av', wspeed_av )
2615	ENDIF
2616
2617	END SUBROUTINE doq_rrd_local_mpi
2618
2619	!--------------------------------------------------------------------------------------------------!
2620	! Description:
2621	! ------------
2622	!> Subroutine writes local (subdomain) restart data
2623	!--------------------------------------------------------------------------------------------------!
2624	SUBROUTINE doq_wrd_local
2625
2626	IMPLICIT NONE
2627
2628
2629	IF ( TRIM( restart_data_format_output ) == 'fortran_binary' ) THEN
2630
2631	IF ( ALLOCATED( pt_2m_av ) ) THEN
2632	CALL wrd_write_string( 'pt_2m_av' )
2633	WRITE ( 14 ) pt_2m_av
2634	ENDIF
2635
2636	IF ( ALLOCATED( rh_av ) ) THEN
2637	CALL wrd_write_string( 'rh_av' )
2638	WRITE ( 14 ) rh_av
2639	ENDIF
2640
2641	IF ( ALLOCATED( ta_av ) ) THEN
2642	CALL wrd_write_string( 'ta_av' )
2643	WRITE ( 14 ) ta_av
2644	ENDIF
2645
2646	IF ( ALLOCATED( ti_av ) ) THEN
2647	CALL wrd_write_string( 'ti_av' )
2648	WRITE ( 14 ) ti_av
2649	ENDIF
2650
2651	IF ( ALLOCATED( u_center_av ) ) THEN
2652	CALL wrd_write_string( 'u_center_av' )
2653	WRITE ( 14 ) u_center_av
2654	ENDIF
2655
2656	IF ( ALLOCATED( uu_av ) ) THEN
2657	CALL wrd_write_string( 'uu_av' )
2658	WRITE ( 14 ) uu_av
2659	ENDIF
2660
2661	IF ( ALLOCATED( uv_10m_av ) ) THEN
2662	CALL wrd_write_string( 'uv_10m_av' )
2663	WRITE ( 14 ) uv_10m_av
2664	ENDIF
2665
2666	IF ( ALLOCATED( v_center_av ) ) THEN
2667	CALL wrd_write_string( 'v_center_av' )
2668	WRITE ( 14 ) v_center_av
2669	ENDIF
2670
2671	IF ( ALLOCATED( vv_av ) ) THEN
2672	CALL wrd_write_string( 'vv_av' )
2673	WRITE ( 14 ) vv_av
2674	ENDIF
2675
2676	IF ( ALLOCATED( ww_av ) ) THEN
2677	CALL wrd_write_string( 'ww_av' )
2678	WRITE ( 14 ) ww_av
2679	ENDIF
2680
2681	IF ( ALLOCATED( wu_av ) ) THEN
2682	CALL wrd_write_string( 'wu_av' )
2683	WRITE ( 14 ) wu_av
2684	ENDIF
2685
2686	IF ( ALLOCATED( wv_av ) ) THEN
2687	CALL wrd_write_string( 'wv_av' )
2688	WRITE ( 14 ) wv_av
2689	ENDIF
2690
2691	IF ( ALLOCATED( wtheta_av ) ) THEN
2692	CALL wrd_write_string( 'wtheta_av' )
2693	WRITE ( 14 ) wtheta_av
2694	ENDIF
2695
2696	IF ( ALLOCATED( wq_av ) ) THEN
2697	CALL wrd_write_string( 'wq_av' )
2698	WRITE ( 14 ) wq_av
2699	ENDIF
2700
2701	IF ( ALLOCATED( ws_av ) ) THEN
2702	CALL wrd_write_string( 'ws_av' )
2703	WRITE ( 14 ) ws_av
2704	ENDIF
2705
2706	IF ( ALLOCATED( wspeed_av ) ) THEN
2707	CALL wrd_write_string( 'wspeed_av' )
2708	WRITE ( 14 ) wspeed_av
2709	ENDIF
2710
2711	ELSEIF ( restart_data_format_output(1:3) == 'mpi' ) THEN
2712
2713	IF ( ALLOCATED( pt_2m_av ) ) CALL wrd_mpi_io( 'pt_2m_av', pt_2m_av )
2714	IF ( ALLOCATED( rh_av ) ) CALL wrd_mpi_io( 'rh_av', rh_av )
2715	IF ( ALLOCATED( ta_av ) ) CALL wrd_mpi_io( 'ta_av', ta_av )
2716	IF ( ALLOCATED( ti_av ) ) CALL wrd_mpi_io( 'ti_av', ti_av )
2717	IF ( ALLOCATED( u_center_av ) ) CALL wrd_mpi_io( 'u_center_av', u_center_av )
2718	IF ( ALLOCATED( uu_av ) ) CALL wrd_mpi_io( 'uu_av', uu_av )
2719	IF ( ALLOCATED( uv_10m_av ) ) CALL wrd_mpi_io( 'uv_10m_av', uv_10m_av )
2720	IF ( ALLOCATED( vv_av ) ) CALL wrd_mpi_io( 'vv_av', vv_av )
2721	IF ( ALLOCATED( v_center_av ) ) CALL wrd_mpi_io( 'v_center_av', v_center_av )
2722	IF ( ALLOCATED( wtheta_av ) ) CALL wrd_mpi_io( 'wtheta_av', wtheta_av )
2723	IF ( ALLOCATED( wq_av ) ) CALL wrd_mpi_io( 'wq_av', wq_av )
2724	IF ( ALLOCATED( ws_av ) ) CALL wrd_mpi_io( 'ws_av', ws_av )
2725	IF ( ALLOCATED( wspeed_av ) ) CALL wrd_mpi_io( 'wspeed_av', wspeed_av )
2726	IF ( ALLOCATED( wu_av ) ) CALL wrd_mpi_io( 'wu_av', wu_av )
2727	IF ( ALLOCATED( wv_av ) ) CALL wrd_mpi_io( 'wv_av', wv_av )
2728	IF ( ALLOCATED( ww_av ) ) CALL wrd_mpi_io( 'ww_av', ww_av )
2729
2730	ENDIF
2731
2732	END SUBROUTINE doq_wrd_local
2733
2734
2735	END MODULE diagnostic_output_quantities_mod

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

Context Navigation

source: palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90

Download in other formats: