Home

Context Navigation

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

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