Home

Context Navigation

source: palm/trunk/SOURCE/diagnostic_output_quantities_mod.f90 @ 4875

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |