Home

Context Navigation

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

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