Home

Context Navigation

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

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