Home

Context Navigation

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

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