Home

Context Navigation

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

Last change on this file since 4861 was 4861, checked in by raasch, 4 years ago
Implemented air temperature as diagnostic output quantity
Property svn:keywords set to `Id`
File size: 93.0 KB

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