Home

Context Navigation

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

Last change on this file since 4759 was 4757, checked in by schwenkel, 4 years ago
implement relative humidity as diagnostic output quantity
Property svn:keywords set to `Id`
File size: 89.3 KB

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