Home

Context Navigation

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

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