Home

Context Navigation

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

Last change on this file since 4458 was 4457, checked in by raasch, 4 years ago
ghost point exchange modularized, bugfix for wrong 2d-exchange
Property svn:keywords set to `Id`
File size: 71.0 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |