Home

Context Navigation

diagnostic_output_quantities_mod.f90 @ 4896

Last change on this file since 4896 was 4895, checked in by suehring, 4 years ago
Remove offset in terrain-following masked output and allow only mask_k_over_surface >= 1
Property svn:keywords set to `Id`
File size: 97.3 KB

Rev	Line
[3994]	1	!> @file diagnostic_output_quantities_mod.f90
[4583]	2	!--------------------------------------------------------------------------------------------------!
[3994]	3	! This file is part of the PALM model system.
	4	!
[4583]	5	! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
	6	! Public License as published by the Free Software Foundation, either version 3 of the License, or
	7	! (at your option) any later version.
[3994]	8	!
[4583]	9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
	10	! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
	11	! Public License for more details.
[3994]	12	!
[4583]	13	! You should have received a copy of the GNU General Public License along with PALM. If not, see
	14	! <http://www.gnu.org/licenses/>.
[3994]	15	!
[4828]	16	! Copyright 1997-2021 Leibniz Universitaet Hannover
[4583]	17	!--------------------------------------------------------------------------------------------------!
[3994]	18	!
	19	! Current revisions:
	20	! ------------------
[4671]	21	!
[4866]	22	!
[3994]	23	! Former revisions:
	24	! -----------------
	25	! $Id: diagnostic_output_quantities_mod.f90 4895 2021-03-03 15:39:08Z raasch $
[4895]	26	! Remove offset in terrain-following masked output
	27	!
	28	! 4866 2021-02-05 17:50:14Z suehring
[4866]	29	! Implemented vertical passive scalar flux
	30	! Implemented air temperature as diagnostic output quantity (previous commit)
	31	!
	32	! 4861 2021-02-01 10:51:45Z raasch
[4768]	33	! Enable 3D data output also with 64-bit precision
	34	!
	35	! 4757 2020-10-26 10:23:38Z schwenkel
[4757]	36	! Implement relative humidity as diagnostic output quantity
	37	!
	38	! 4671 2020-09-09 20:27:58Z pavelkrc
[4671]	39	! Implementation of downward facing USM and LSM surfaces
[4861]	40	!
[4671]	41	! 4583 2020-06-29 12:36:47Z raasch
[4583]	42	! file re-formatted to follow the PALM coding standard
	43	!
	44	! 4560 2020-06-11 12:19:47Z suehring
[4560]	45	! - Bugfix in calculation of vertical momentum and scalar fluxes
	46	! - remove averaged output variables from PUBLIC list
[4583]	47	!
[4560]	48	! 4535 2020-05-15 12:07:23Z raasch
[4535]	49	! bugfix for restart data format query
[4583]	50	!
[4535]	51	! 4518 2020-05-04 15:44:28Z suehring
[4583]	52	! * Define arrays over ghost points in order to allow for standard mpi-io treatment. By this
	53	! modularization of restart-data input is possible with the module interface.
[4518]	54	! * Bugfix: add missing restart input of wtheta_av, wq_av, wu_av, and wv_av.
[4583]	55	!
[4518]	56	! 4517 2020-05-03 14:29:30Z raasch
[4457]	57	! use statement for exchange horiz added,
	58	! bugfix for call of exchange horiz 2d
[4583]	59	!
[4457]	60	! 4431 2020-02-27 23:23:01Z gronemeier
[4431]	61	! added wspeed and wdir output; bugfix: set fill_value in case of masked output
	62	!
	63	! 4360 2020-01-07 11:25:50Z suehring
[4583]	64	! added output of wu, wv, wtheta and wq to enable covariance calculation according to temporal EC
	65	! method
[4431]	66	!
[4351]	67	! 4346 2019-12-18 11:55:56Z motisi
[4583]	68	! Introduction of wall_flags_total_0, which currently sets bits based on static topography
	69	! information used in wall_flags_static_0
[4431]	70	!
[4346]	71	! 4331 2019-12-10 18:25:02Z suehring
[4331]	72	! - Modularize 2-m potential temperature output
	73	! - New output for 10-m wind speed
[4431]	74	!
[4331]	75	! 4329 2019-12-10 15:46:36Z motisi
[4329]	76	! Renamed wall_flags_0 to wall_flags_static_0
[4431]	77	!
[4329]	78	! 4182 2019-08-22 15:20:23Z scharf
[4182]	79	! Corrected "Former revisions" section
[4431]	80	!
[4182]	81	! 4167 2019-08-16 11:01:48Z suehring
[4583]	82	! Changed behaviour of masked output over surface to follow terrain and ignore buildings
	83	! (J.Resler, T.Gronemeier)
[4431]	84	!
[4167]	85	! 4157 2019-08-14 09:19:12Z suehring
[4583]	86	! Initialization restructured, in order to work also when data output during spin-up is enabled.
[4431]	87	!
[4157]	88	! 4132 2019-08-02 12:34:17Z suehring
[4132]	89	! Bugfix in masked data output
[4431]	90	!
[4132]	91	! 4069 2019-07-01 14:05:51Z Giersch
[4583]	92	! Masked output running index mid has been introduced as a local variable to avoid runtime error
	93	! (Loop variable has been modified) in time_integration
[4431]	94	!
[4069]	95	! 4039 2019-06-18 10:32:41Z suehring
[4583]	96	! - Add output of uu, vv, ww to enable variance calculation according temporal EC method
[4039]	97	! - Allocate arrays only when they are required
	98	! - Formatting adjustment
	99	! - Rename subroutines
	100	! - Further modularization
[4431]	101	!
[4039]	102	! 3998 2019-05-23 13:38:11Z suehring
[4431]	103	! Bugfix in gathering all output strings
	104	!
[3998]	105	! 3995 2019-05-22 18:59:54Z suehring
[4583]	106	! Avoid compiler warnings about unused variable and fix string operation which is not allowed with
	107	! PGI compiler
[4431]	108	!
[3995]	109	! 3994 2019-05-22 18:08:09Z suehring
[3994]	110	! Initial revision
	111	!
[4182]	112	! Authors:
	113	! --------
[3994]	114	! @author Farah Kanani-Suehring
[4182]	115	!
[4431]	116	!
[3994]	117	! Description:
	118	! ------------
	119	!> ...
[4583]	120	!--------------------------------------------------------------------------------------------------!
[3994]	121	MODULE diagnostic_output_quantities_mod
[4431]	122
[4583]	123	USE arrays_3d, &
	124	ONLY: ddzu, &
[4757]	125	exner, &
	126	hyp, &
[4583]	127	pt, &
	128	q, &
[4866]	129	s, &
[4583]	130	u, &
	131	v, &
	132	w, &
	133	zu, &
[4331]	134	zw
	135
[4583]	136	USE basic_constants_and_equations_mod, &
[4861]	137	ONLY: degc_to_k, &
	138	kappa, &
[4757]	139	pi, &
	140	magnus, &
	141	rd_d_rv
[4331]	142
[4583]	143	USE control_parameters, &
	144	ONLY: current_timestep_number, &
	145	data_output, &
[4757]	146	data_output_pr, &
	147	humidity, &
[4583]	148	message_string, &
[4866]	149	passive_scalar, &
[4583]	150	restart_data_format_output, &
[4331]	151	varnamelength
[4431]	152	!
[4583]	153	! USE cpulog, &
[3994]	154	! ONLY: cpu_log, log_point
[4039]	155
[4583]	156	USE exchange_horiz_mod, &
[4457]	157	ONLY: exchange_horiz_2d
	158
[4583]	159	USE grid_variables, &
[4039]	160	ONLY: ddx, ddy
[4331]	161
[4583]	162	USE indices, &
	163	ONLY: nbgp, &
	164	nxl, &
	165	nxlg, &
	166	nxr, &
	167	nxrg, &
	168	nyn, &
	169	nyng, &
	170	nys, &
	171	nysg, &
	172	nzb, &
	173	nzt, &
[4346]	174	wall_flags_total_0
[4331]	175
[3994]	176	USE kinds
	177
[4757]	178	USE profil_parameter, &
	179	ONLY: dopr_index
	180
[4518]	181	USE restart_data_mpi_io_mod, &
	182	ONLY: rd_mpi_io_check_array, &
	183	rrd_mpi_io, &
	184	wrd_mpi_io
	185
[4757]	186	USE statistics, &
	187	ONLY: hom, &
	188	pr_palm, &
	189	rmask, &
	190	statistic_regions, &
	191	sums_l
	192
[4583]	193	USE surface_mod, &
	194	ONLY: surf_def_h, &
	195	surf_lsm_h, &
	196	surf_type, &
[4331]	197	surf_usm_h
[3994]	198
	199	IMPLICIT NONE
	200
	201	CHARACTER(LEN=varnamelength), DIMENSION(500) :: do_all = ' '
	202
	203	INTEGER(iwp) :: timestep_number_at_prev_calc = 0 !< ...at previous diagnostic output calculation
[4431]	204
[4039]	205	LOGICAL :: initialized_diagnostic_output_quantities = .FALSE. !< flag indicating whether output is initialized
	206	LOGICAL :: prepared_diagnostic_output_quantities = .FALSE. !< flag indicating whether output is p
[3994]	207
[4331]	208	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pt_2m !< 2-m air potential temperature
	209	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: pt_2m_av !< averaged 2-m air potential temperature
	210	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: uv_10m !< horizontal wind speed at 10m
	211	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: uv_10m_av !< averaged horizontal wind speed at 10m
	212
[4757]	213	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: rh !< relative humidity
	214	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: rh_av !< avg. relative humidity
[4861]	215	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ta !< air temperature
	216	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ta_av !< avg. air temperature
[3994]	217	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ti !< rotation(u,v,w) aka turbulence intensity
[4431]	218	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ti_av !< avg. rotation(u,v,w) aka turbulence intensity
	219	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: u_center !< u at center of grid box
	220	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: u_center_av !< mean of u_center
	221	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: uu !< uu
	222	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: uu_av !< mean of uu
	223	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: wspeed !< horizontal wind speed
	224	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: wspeed_av !< mean of horizotal wind speed
	225	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: v_center !< v at center of grid box
	226	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: v_center_av !< mean of v_center
	227	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: vv !< vv
	228	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: vv_av !< mean of vv
	229	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: wdir !< wind direction
	230	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: wdir_av !< mean wind direction
	231	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ww !< ww
	232	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ww_av !< mean of ww
	233	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: wu !< wu
	234	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: wu_av !< mean of wu
	235	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: wv !< wv
	236	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: wv_av !< mean of wv
	237	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: wtheta !< wtheta
	238	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: wtheta_av !< mean of wtheta
	239	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: wq !< wq
	240	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: wq_av !< mean of wq
[4866]	241	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ws !< ws
	242	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ws_av !< mean of ws
[3994]	243
	244
	245	SAVE
	246
	247	PRIVATE
	248
	249	!
	250	!-- Public variables
[4583]	251	PUBLIC do_all, &
	252	initialized_diagnostic_output_quantities, &
	253	prepared_diagnostic_output_quantities, &
[4560]	254	timestep_number_at_prev_calc
[4431]	255	!
	256	!-- Public routines
[4583]	257	PUBLIC doq_3d_data_averaging, &
	258	doq_calculate, &
	259	doq_check_data_output, &
[4757]	260	doq_check_data_output_pr, &
[4583]	261	doq_define_netcdf_grid, &
	262	doq_init, &
	263	doq_output_2d, &
	264	doq_output_3d, &
	265	doq_output_mask, &
[4757]	266	doq_statistics, &
[4583]	267	doq_rrd_local, &
[4039]	268	doq_wrd_local
[3994]	269
	270
[4039]	271	INTERFACE doq_3d_data_averaging
	272	MODULE PROCEDURE doq_3d_data_averaging
[4431]	273	END INTERFACE doq_3d_data_averaging
[3994]	274
[4039]	275	INTERFACE doq_calculate
	276	MODULE PROCEDURE doq_calculate
	277	END INTERFACE doq_calculate
[3994]	278
[4039]	279	INTERFACE doq_check_data_output
	280	MODULE PROCEDURE doq_check_data_output
	281	END INTERFACE doq_check_data_output
[4431]	282
[4757]	283	INTERFACE doq_check_data_output_pr
	284	MODULE PROCEDURE doq_check_data_output_pr
	285	END INTERFACE doq_check_data_output_pr
	286
[4039]	287	INTERFACE doq_define_netcdf_grid
	288	MODULE PROCEDURE doq_define_netcdf_grid
	289	END INTERFACE doq_define_netcdf_grid
[4431]	290
[4039]	291	INTERFACE doq_output_2d
	292	MODULE PROCEDURE doq_output_2d
	293	END INTERFACE doq_output_2d
[4431]	294
[4039]	295	INTERFACE doq_output_3d
	296	MODULE PROCEDURE doq_output_3d
	297	END INTERFACE doq_output_3d
[4431]	298
[4039]	299	INTERFACE doq_output_mask
	300	MODULE PROCEDURE doq_output_mask
	301	END INTERFACE doq_output_mask
[4431]	302
[4039]	303	INTERFACE doq_init
	304	MODULE PROCEDURE doq_init
	305	END INTERFACE doq_init
[3994]	306
[4757]	307	INTERFACE doq_statistics
	308	MODULE PROCEDURE doq_statistics
	309	END INTERFACE doq_statistics
	310
[4039]	311	INTERFACE doq_prepare
	312	MODULE PROCEDURE doq_prepare
	313	END INTERFACE doq_prepare
[4431]	314
[4518]	315	INTERFACE doq_rrd_local
	316	MODULE PROCEDURE doq_rrd_local_ftn
	317	MODULE PROCEDURE doq_rrd_local_mpi
	318	END INTERFACE doq_rrd_local
[4431]	319
[4039]	320	INTERFACE doq_wrd_local
	321	MODULE PROCEDURE doq_wrd_local
	322	END INTERFACE doq_wrd_local
[3994]	323
	324	CONTAINS
[4431]	325
[4583]	326	!--------------------------------------------------------------------------------------------------!
[4039]	327	! Description:
	328	! ------------
[4583]	329	!> Sum up and time-average diagnostic output quantities as well as allocate the array necessary for
	330	!> storing the average.
	331	!--------------------------------------------------------------------------------------------------!
[4039]	332	SUBROUTINE doq_3d_data_averaging( mode, variable )
[3994]	333
[4583]	334	USE control_parameters, &
[4039]	335	ONLY: average_count_3d
	336
[4431]	337	CHARACTER (LEN=*) :: mode !<
	338	CHARACTER (LEN=*) :: variable !<
[4039]	339
	340	INTEGER(iwp) :: i !<
	341	INTEGER(iwp) :: j !<
	342	INTEGER(iwp) :: k !<
	343
[4583]	344
[4039]	345	IF ( mode == 'allocate' ) THEN
	346
	347	SELECT CASE ( TRIM( variable ) )
	348
[4757]	349	CASE ( 'rh' )
	350	IF ( .NOT. ALLOCATED( rh_av ) ) THEN
	351	ALLOCATE( rh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	352	ENDIF
	353	rh_av = 0.0_wp
	354
[4861]	355	CASE ( 'ta' )
	356	IF ( .NOT. ALLOCATED( ta_av ) ) THEN
	357	ALLOCATE( ta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	358	ENDIF
	359	ta_av = 0.0_wp
	360
[4039]	361	CASE ( 'ti' )
	362	IF ( .NOT. ALLOCATED( ti_av ) ) THEN
[4518]	363	ALLOCATE( ti_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4039]	364	ENDIF
	365	ti_av = 0.0_wp
[4431]	366
[4039]	367	CASE ( 'uu' )
	368	IF ( .NOT. ALLOCATED( uu_av ) ) THEN
[4518]	369	ALLOCATE( uu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4039]	370	ENDIF
	371	uu_av = 0.0_wp
[4431]	372
[4039]	373	CASE ( 'vv' )
	374	IF ( .NOT. ALLOCATED( vv_av ) ) THEN
[4518]	375	ALLOCATE( vv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4039]	376	ENDIF
	377	vv_av = 0.0_wp
[4431]	378
[4039]	379	CASE ( 'ww' )
	380	IF ( .NOT. ALLOCATED( ww_av ) ) THEN
[4518]	381	ALLOCATE( ww_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4039]	382	ENDIF
	383	ww_av = 0.0_wp
[4431]	384
[4351]	385	CASE ( 'wu' )
	386	IF ( .NOT. ALLOCATED( wu_av ) ) THEN
[4518]	387	ALLOCATE( wu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4351]	388	ENDIF
	389	wu_av = 0.0_wp
[4431]	390
[4351]	391	CASE ( 'wv' )
	392	IF ( .NOT. ALLOCATED( wv_av ) ) THEN
[4518]	393	ALLOCATE( wv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4351]	394	ENDIF
	395	wv_av = 0.0_wp
[4431]	396
[4351]	397	CASE ( 'wtheta' )
	398	IF ( .NOT. ALLOCATED( wtheta_av ) ) THEN
[4518]	399	ALLOCATE( wtheta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4351]	400	ENDIF
	401	wtheta_av = 0.0_wp
[4431]	402
[4351]	403	CASE ( 'wq' )
	404	IF ( .NOT. ALLOCATED( wq_av ) ) THEN
[4518]	405	ALLOCATE( wq_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4351]	406	ENDIF
	407	wq_av = 0.0_wp
[4431]	408
[4866]	409	CASE ( 'ws' )
	410	IF ( .NOT. ALLOCATED( ws_av ) ) THEN
	411	ALLOCATE( ws_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	412	ENDIF
	413	ws_av = 0.0_wp
	414
[4331]	415	CASE ( 'theta_2m*' )
	416	IF ( .NOT. ALLOCATED( pt_2m_av ) ) THEN
	417	ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )
	418	ENDIF
	419	pt_2m_av = 0.0_wp
[4431]	420
[4331]	421	CASE ( 'wspeed_10m*' )
	422	IF ( .NOT. ALLOCATED( uv_10m_av ) ) THEN
	423	ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) )
	424	ENDIF
	425	uv_10m_av = 0.0_wp
[4039]	426
[4431]	427	CASE ( 'wspeed' )
	428	IF ( .NOT. ALLOCATED( wspeed_av ) ) THEN
[4518]	429	ALLOCATE( wspeed_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4431]	430	ENDIF
	431	wspeed_av = 0.0_wp
	432
	433	CASE ( 'wdir' )
	434	IF ( .NOT. ALLOCATED( u_center_av ) ) THEN
[4518]	435	ALLOCATE( u_center_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4431]	436	ENDIF
	437	IF ( .NOT. ALLOCATED( v_center_av ) ) THEN
[4518]	438	ALLOCATE( v_center_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4431]	439	ENDIF
	440	u_center_av = 0.0_wp
	441	v_center_av = 0.0_wp
	442
[4039]	443	CASE DEFAULT
	444	CONTINUE
	445
	446	END SELECT
	447
	448	ELSEIF ( mode == 'sum' ) THEN
	449
	450	SELECT CASE ( TRIM( variable ) )
[4431]	451
[4757]	452	CASE ( 'rh' )
	453	IF ( ALLOCATED( rh_av ) ) THEN
	454	DO i = nxl, nxr
	455	DO j = nys, nyn
	456	DO k = nzb, nzt+1
	457	rh_av(k,j,i) = rh_av(k,j,i) + rh(k,j,i)
	458	ENDDO
	459	ENDDO
	460	ENDDO
	461	ENDIF
	462
[4861]	463	CASE ( 'ta' )
	464	IF ( ALLOCATED( ta_av ) ) THEN
	465	DO i = nxl, nxr
	466	DO j = nys, nyn
	467	DO k = nzb, nzt+1
	468	ta_av(k,j,i) = ta_av(k,j,i) + ta(k,j,i)
	469	ENDDO
	470	ENDDO
	471	ENDDO
	472	ENDIF
	473
[4039]	474	CASE ( 'ti' )
[4583]	475	IF ( ALLOCATED( ti_av ) ) THEN
[4039]	476	DO i = nxl, nxr
	477	DO j = nys, nyn
	478	DO k = nzb, nzt+1
	479	ti_av(k,j,i) = ti_av(k,j,i) + ti(k,j,i)
	480	ENDDO
	481	ENDDO
	482	ENDDO
	483	ENDIF
[4431]	484
[4039]	485	CASE ( 'uu' )
[4583]	486	IF ( ALLOCATED( uu_av ) ) THEN
[4039]	487	DO i = nxl, nxr
	488	DO j = nys, nyn
	489	DO k = nzb, nzt+1
	490	uu_av(k,j,i) = uu_av(k,j,i) + uu(k,j,i)
	491	ENDDO
	492	ENDDO
	493	ENDDO
	494	ENDIF
[4431]	495
[4039]	496	CASE ( 'vv' )
[4583]	497	IF ( ALLOCATED( vv_av ) ) THEN
[4039]	498	DO i = nxl, nxr
	499	DO j = nys, nyn
	500	DO k = nzb, nzt+1
	501	vv_av(k,j,i) = vv_av(k,j,i) + vv(k,j,i)
	502	ENDDO
	503	ENDDO
	504	ENDDO
	505	ENDIF
[4431]	506
[4039]	507	CASE ( 'ww' )
[4583]	508	IF ( ALLOCATED( ww_av ) ) THEN
[4039]	509	DO i = nxl, nxr
	510	DO j = nys, nyn
	511	DO k = nzb, nzt+1
	512	ww_av(k,j,i) = ww_av(k,j,i) + ww(k,j,i)
	513	ENDDO
	514	ENDDO
	515	ENDDO
	516	ENDIF
[4431]	517
[4351]	518	CASE ( 'wu' )
[4583]	519	IF ( ALLOCATED( wu_av ) ) THEN
[4351]	520	DO i = nxl, nxr
	521	DO j = nys, nyn
	522	DO k = nzb, nzt+1
	523	wu_av(k,j,i) = wu_av(k,j,i) + wu(k,j,i)
	524	ENDDO
	525	ENDDO
	526	ENDDO
	527	ENDIF
[4431]	528
[4351]	529	CASE ( 'wv' )
[4583]	530	IF ( ALLOCATED( wv_av ) ) THEN
[4351]	531	DO i = nxl, nxr
	532	DO j = nys, nyn
	533	DO k = nzb, nzt+1
	534	wv_av(k,j,i) = wv_av(k,j,i) + wv(k,j,i)
	535	ENDDO
	536	ENDDO
	537	ENDDO
	538	ENDIF
[4431]	539
[4351]	540	CASE ( 'wtheta' )
[4583]	541	IF ( ALLOCATED( wtheta_av ) ) THEN
[4351]	542	DO i = nxl, nxr
	543	DO j = nys, nyn
	544	DO k = nzb, nzt+1
	545	wtheta_av(k,j,i) = wtheta_av(k,j,i) + wtheta(k,j,i)
	546	ENDDO
	547	ENDDO
	548	ENDDO
	549	ENDIF
[4431]	550
[4351]	551	CASE ( 'wq' )
[4583]	552	IF ( ALLOCATED( wq_av ) ) THEN
[4351]	553	DO i = nxl, nxr
	554	DO j = nys, nyn
	555	DO k = nzb, nzt+1
	556	wq_av(k,j,i) = wq_av(k,j,i) + wq(k,j,i)
	557	ENDDO
	558	ENDDO
	559	ENDDO
	560	ENDIF
[4431]	561
[4866]	562	CASE ( 'ws' )
	563	IF ( ALLOCATED( ws_av ) ) THEN
	564	DO i = nxl, nxr
	565	DO j = nys, nyn
	566	DO k = nzb, nzt+1
	567	ws_av(k,j,i) = ws_av(k,j,i) + ws(k,j,i)
	568	ENDDO
	569	ENDDO
	570	ENDDO
	571	ENDIF
	572
[4331]	573	CASE ( 'theta_2m*' )
[4583]	574	IF ( ALLOCATED( pt_2m_av ) ) THEN
[4331]	575	DO i = nxl, nxr
	576	DO j = nys, nyn
	577	pt_2m_av(j,i) = pt_2m_av(j,i) + pt_2m(j,i)
	578	ENDDO
	579	ENDDO
	580	ENDIF
[4039]	581
[4331]	582	CASE ( 'wspeed_10m*' )
[4583]	583	IF ( ALLOCATED( uv_10m_av ) ) THEN
[4331]	584	DO i = nxl, nxr
	585	DO j = nys, nyn
	586	uv_10m_av(j,i) = uv_10m_av(j,i) + uv_10m(j,i)
	587	ENDDO
	588	ENDDO
	589	ENDIF
	590
[4431]	591	CASE ( 'wspeed' )
[4583]	592	IF ( ALLOCATED( wspeed_av ) ) THEN
[4431]	593	DO i = nxl, nxr
	594	DO j = nys, nyn
	595	DO k = nzb, nzt+1
	596	wspeed_av(k,j,i) = wspeed_av(k,j,i) + wspeed(k,j,i)
	597	ENDDO
	598	ENDDO
	599	ENDDO
	600	ENDIF
	601
	602	CASE ( 'wdir' )
[4583]	603	IF ( ALLOCATED( u_center_av ) .AND. ALLOCATED( v_center_av ) ) THEN
[4431]	604	DO i = nxl, nxr
	605	DO j = nys, nyn
	606	DO k = nzb, nzt+1
	607	u_center_av(k,j,i) = u_center_av(k,j,i) + u_center(k,j,i)
	608	v_center_av(k,j,i) = v_center_av(k,j,i) + v_center(k,j,i)
	609	ENDDO
	610	ENDDO
	611	ENDDO
	612	ENDIF
	613
[4039]	614	CASE DEFAULT
	615	CONTINUE
	616
	617	END SELECT
	618
	619	ELSEIF ( mode == 'average' ) THEN
	620
	621	SELECT CASE ( TRIM( variable ) )
	622
[4757]	623	CASE ( 'rh' )
	624	IF ( ALLOCATED( rh_av ) ) THEN
	625	DO i = nxl, nxr
	626	DO j = nys, nyn
	627	DO k = nzb, nzt+1
	628	rh_av(k,j,i) = rh_av(k,j,i) / REAL( average_count_3d, KIND=wp )
	629	ENDDO
	630	ENDDO
	631	ENDDO
	632	ENDIF
	633
[4861]	634	CASE ( 'ta' )
	635	IF ( ALLOCATED( ta_av ) ) THEN
	636	DO i = nxl, nxr
	637	DO j = nys, nyn
	638	DO k = nzb, nzt+1
	639	ta_av(k,j,i) = ta_av(k,j,i) / REAL( average_count_3d, KIND=wp )
	640	ENDDO
	641	ENDDO
	642	ENDDO
	643	ENDIF
	644
[4039]	645	CASE ( 'ti' )
[4583]	646	IF ( ALLOCATED( ti_av ) ) THEN
[4039]	647	DO i = nxl, nxr
	648	DO j = nys, nyn
	649	DO k = nzb, nzt+1
	650	ti_av(k,j,i) = ti_av(k,j,i) / REAL( average_count_3d, KIND=wp )
	651	ENDDO
	652	ENDDO
	653	ENDDO
	654	ENDIF
[4431]	655
[4039]	656	CASE ( 'uu' )
[4583]	657	IF ( ALLOCATED( uu_av ) ) THEN
[4039]	658	DO i = nxl, nxr
	659	DO j = nys, nyn
	660	DO k = nzb, nzt+1
	661	uu_av(k,j,i) = uu_av(k,j,i) / REAL( average_count_3d, KIND=wp )
	662	ENDDO
	663	ENDDO
	664	ENDDO
	665	ENDIF
[4431]	666
[4039]	667	CASE ( 'vv' )
[4583]	668	IF ( ALLOCATED( vv_av ) ) THEN
[4039]	669	DO i = nxl, nxr
	670	DO j = nys, nyn
	671	DO k = nzb, nzt+1
	672	vv_av(k,j,i) = vv_av(k,j,i) / REAL( average_count_3d, KIND=wp )
	673	ENDDO
	674	ENDDO
	675	ENDDO
	676	ENDIF
[4431]	677
[4039]	678	CASE ( 'ww' )
[4583]	679	IF ( ALLOCATED( ww_av ) ) THEN
[4039]	680	DO i = nxl, nxr
	681	DO j = nys, nyn
	682	DO k = nzb, nzt+1
	683	ww_av(k,j,i) = ww_av(k,j,i) / REAL( average_count_3d, KIND=wp )
	684	ENDDO
	685	ENDDO
	686	ENDDO
	687	ENDIF
	688
[4351]	689	CASE ( 'wu' )
[4583]	690	IF ( ALLOCATED( wu_av ) ) THEN
[4351]	691	DO i = nxl, nxr
	692	DO j = nys, nyn
	693	DO k = nzb, nzt+1
	694	wu_av(k,j,i) = wu_av(k,j,i) / REAL( average_count_3d, KIND=wp )
	695	ENDDO
	696	ENDDO
	697	ENDDO
	698	ENDIF
[4431]	699
[4351]	700	CASE ( 'wv' )
[4583]	701	IF ( ALLOCATED( wv_av ) ) THEN
[4351]	702	DO i = nxl, nxr
	703	DO j = nys, nyn
	704	DO k = nzb, nzt+1
	705	wv_av(k,j,i) = wv_av(k,j,i) / REAL( average_count_3d, KIND=wp )
	706	ENDDO
	707	ENDDO
	708	ENDDO
	709	ENDIF
[4431]	710
[4351]	711	CASE ( 'wtheta' )
[4583]	712	IF ( ALLOCATED( wtheta_av ) ) THEN
[4351]	713	DO i = nxl, nxr
	714	DO j = nys, nyn
	715	DO k = nzb, nzt+1
	716	wtheta_av(k,j,i) = wtheta_av(k,j,i) / REAL( average_count_3d, KIND=wp )
	717	ENDDO
	718	ENDDO
	719	ENDDO
	720	ENDIF
[4431]	721
[4351]	722	CASE ( 'wq' )
[4583]	723	IF ( ALLOCATED( wq_av ) ) THEN
[4351]	724	DO i = nxl, nxr
	725	DO j = nys, nyn
	726	DO k = nzb, nzt+1
	727	wq_av(k,j,i) = wq_av(k,j,i) / REAL( average_count_3d, KIND=wp )
	728	ENDDO
	729	ENDDO
	730	ENDDO
	731	ENDIF
	732
[4866]	733	CASE ( 'ws' )
	734	IF ( ALLOCATED( ws_av ) ) THEN
	735	DO i = nxl, nxr
	736	DO j = nys, nyn
	737	DO k = nzb, nzt+1
	738	ws_av(k,j,i) = ws_av(k,j,i) / REAL( average_count_3d, KIND=wp )
	739	ENDDO
	740	ENDDO
	741	ENDDO
	742	ENDIF
	743
[4331]	744	CASE ( 'theta_2m*' )
[4583]	745	IF ( ALLOCATED( pt_2m_av ) ) THEN
[4331]	746	DO i = nxlg, nxrg
	747	DO j = nysg, nyng
	748	pt_2m_av(j,i) = pt_2m_av(j,i) / REAL( average_count_3d, KIND=wp )
	749	ENDDO
	750	ENDDO
[4457]	751	CALL exchange_horiz_2d( pt_2m_av )
[4331]	752	ENDIF
	753
	754	CASE ( 'wspeed_10m*' )
[4583]	755	IF ( ALLOCATED( uv_10m_av ) ) THEN
[4331]	756	DO i = nxlg, nxrg
	757	DO j = nysg, nyng
	758	uv_10m_av(j,i) = uv_10m_av(j,i) / REAL( average_count_3d, KIND=wp )
	759	ENDDO
	760	ENDDO
[4457]	761	CALL exchange_horiz_2d( uv_10m_av )
[4331]	762	ENDIF
	763
[4431]	764	CASE ( 'wspeed' )
[4583]	765	IF ( ALLOCATED( wspeed_av ) ) THEN
[4431]	766	DO i = nxl, nxr
	767	DO j = nys, nyn
	768	DO k = nzb, nzt+1
	769	wspeed_av(k,j,i) = wspeed_av(k,j,i) / REAL( average_count_3d, KIND=wp )
	770	ENDDO
	771	ENDDO
	772	ENDDO
	773	ENDIF
	774
	775	CASE ( 'wdir' )
[4583]	776	IF ( ALLOCATED( u_center_av ) .AND. ALLOCATED( v_center_av ) ) THEN
[4431]	777
	778	IF ( .NOT. ALLOCATED( wdir_av ) ) THEN
[4518]	779	ALLOCATE( wdir_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4431]	780	ENDIF
	781	wdir_av = 0.0_wp
	782
	783	DO i = nxl, nxr
	784	DO j = nys, nyn
	785	DO k = nzb, nzt+1
	786	u_center_av(k,j,i) = u_center_av(k,j,i) / REAL( average_count_3d, KIND=wp )
	787	v_center_av(k,j,i) = v_center_av(k,j,i) / REAL( average_count_3d, KIND=wp )
[4583]	788	wdir_av(k,j,i) = ATAN2( u_center_av(k,j,i), v_center_av(k,j,i) ) &
	789	/ pi * 180.0_wp + 180.0_wp
[4431]	790	ENDDO
	791	ENDDO
	792	ENDDO
	793	ENDIF
	794
[4039]	795	END SELECT
	796
	797	ENDIF
	798
	799
[4431]	800	END SUBROUTINE doq_3d_data_averaging
	801
[4583]	802	!--------------------------------------------------------------------------------------------------!
[3994]	803	! Description:
	804	! ------------
[4039]	805	!> Check data output for diagnostic output
[4583]	806	!--------------------------------------------------------------------------------------------------!
[4331]	807	SUBROUTINE doq_check_data_output( var, unit, i, ilen, k )
[3994]	808
[4039]	809	IMPLICIT NONE
[3994]	810
[4431]	811	CHARACTER (LEN=*) :: unit !<
[4039]	812	CHARACTER (LEN=*) :: var !<
[3994]	813
[4331]	814	INTEGER(iwp), OPTIONAL, INTENT(IN) :: i !< Current element of data_output
	815	INTEGER(iwp), OPTIONAL, INTENT(IN) :: ilen !< Length of current entry in data_output
	816	INTEGER(iwp), OPTIONAL, INTENT(IN) :: k !< Output is xy mode? 0 = no, 1 = yes
	817
[4583]	818
[4039]	819	SELECT CASE ( TRIM( var ) )
	820
[4757]	821	CASE ( 'rh' )
	822	unit = '%'
	823
[4861]	824	CASE ( 'ta' )
	825	unit = 'degree_C'
	826
[4039]	827	CASE ( 'ti' )
	828	unit = '1/s'
[4431]	829
[4039]	830	CASE ( 'uu' )
	831	unit = 'm2/s2'
[4431]	832
[4039]	833	CASE ( 'vv' )
	834	unit = 'm2/s2'
[4431]	835
[4039]	836	CASE ( 'ww' )
	837	unit = 'm2/s2'
[4431]	838
[4351]	839	CASE ( 'wu' )
	840	unit = 'm2/s2'
[4431]	841
[4351]	842	CASE ( 'wv' )
	843	unit = 'm2/s2'
[4431]	844
[4351]	845	CASE ( 'wtheta' )
	846	unit = 'Km/s'
[4431]	847
[4351]	848	CASE ( 'wq' )
	849	unit = 'm/s'
[4431]	850
[4866]	851	CASE ( 'ws' )
	852	IF ( .NOT. passive_scalar ) THEN
	853	message_string = 'data_output = ' // TRIM( var ) // &
	854	' is not implemented for passive_scalar = .FALSE.'
	855	CALL message( 'diagnostic_output', 'PA0093', 1, 2, 0, 6, 0 )
	856	ENDIF
	857	unit = 'm/s'
	858
[4431]	859	CASE ( 'wspeed' )
	860	unit = 'm/s'
	861
	862	CASE ( 'wdir' )
	863	unit = 'degree'
[4331]	864	!
	865	!-- Treat horizotal cross-section output quanatities
	866	CASE ( 'theta_2m', 'wspeed_10m' )
	867	!
	868	!-- Check if output quantity is _xy only.
	869	IF ( k == 0 .OR. data_output(i)(ilen-2:ilen) /= '_xy' ) THEN
[4583]	870	message_string = 'illegal value for data_output: "' // &
	871	TRIM( var ) // '" & only 2d-horizontal ' // &
[4331]	872	'cross sections are allowed for this value'
	873	CALL message( 'diagnostic_output', 'PA0111', 1, 2, 0, 6, 0 )
	874	ENDIF
	875
	876	IF ( TRIM( var ) == 'theta_2m*' ) unit = 'K'
	877	IF ( TRIM( var ) == 'wspeed_10m*' ) unit = 'm/s'
	878
[4039]	879	CASE DEFAULT
	880	unit = 'illegal'
	881
	882	END SELECT
	883
	884
	885	END SUBROUTINE doq_check_data_output
[4431]	886
[4757]	887
[4583]	888	!--------------------------------------------------------------------------------------------------!
[4757]	889	! Description:
	890	! ------------
	891	!> Set the unit of user defined profile output quantities. For those variables not recognized by the
	892	!> user, the parameter unit is set to "illegal", which tells the calling routine that the
	893	!> output variable is not defined and leads to a program abort.
	894	!--------------------------------------------------------------------------------------------------!
	895	SUBROUTINE doq_check_data_output_pr( variable, var_count, unit, dopr_unit )
	896
	897	CHARACTER (LEN=*) :: unit !<
	898	CHARACTER (LEN=*) :: variable !<
	899	CHARACTER (LEN=*) :: dopr_unit !< local value of dopr_unit
	900
	901	INTEGER(iwp) :: pr_index !<
	902	INTEGER(iwp) :: var_count !<
	903
	904	SELECT CASE ( TRIM( variable ) )
	905
	906	CASE ( 'rh' )
	907	IF ( .NOT. humidity ) THEN
	908	message_string = 'data_output_pr = ' // TRIM( data_output_pr(var_count) ) // &
	909	' requires humidity'
	910	CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
	911	ENDIF
	912	pr_index = 130
	913	dopr_index(var_count) = pr_index
	914	dopr_unit = '%'
	915	unit = dopr_unit
	916	hom(:,2,pr_index,:) = SPREAD( zu, 2, statistic_regions+1 )
	917
	918	CASE DEFAULT
	919	unit = 'illegal'
	920
	921	END SELECT
	922
	923
	924	END SUBROUTINE doq_check_data_output_pr
	925
	926
	927	!--------------------------------------------------------------------------------------------------!
[4039]	928	!
	929	! Description:
	930	! ------------
	931	!> Subroutine defining appropriate grid for netcdf variables.
[4583]	932	!--------------------------------------------------------------------------------------------------!
[4039]	933	SUBROUTINE doq_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z )
[4431]	934
[3994]	935	IMPLICIT NONE
[4039]	936
	937	CHARACTER (LEN=*), INTENT(IN) :: variable !<
[4431]	938	CHARACTER (LEN=*), INTENT(OUT) :: grid_x !<
	939	CHARACTER (LEN=*), INTENT(OUT) :: grid_y !<
	940	CHARACTER (LEN=*), INTENT(OUT) :: grid_z !<
[4039]	941
[4583]	942	LOGICAL, INTENT(OUT) :: found !<
	943
	944
[4039]	945	found = .TRUE.
[4431]	946
[4039]	947	SELECT CASE ( TRIM( variable ) )
[3995]	948	!
[4039]	949	!-- s grid
[4757]	950	CASE ( 'rh', 'rh_xy', 'rh_xz', 'rh_yz', &
[4861]	951	'ta', 'ta_xy', 'ta_xz', 'ta_yz', &
[4757]	952	'ti', 'ti_xy', 'ti_xz', 'ti_yz', &
[4583]	953	'wspeed', 'wspeed_xy', 'wspeed_xz', 'wspeed_yz', &
[4560]	954	'wdir', 'wdir_xy', 'wdir_xz', 'wdir_yz' )
[3994]	955
[4039]	956	grid_x = 'x'
	957	grid_y = 'y'
[4331]	958	grid_z = 'zu'
[4039]	959	!
[4331]	960	!-- s grid surface variables
	961	CASE ( 'theta_2m_xy', 'wspeed_10m' )
	962
	963	grid_x = 'x'
	964	grid_y = 'y'
	965	grid_z = 'zu'
	966	!
[4039]	967	!-- u grid
	968	CASE ( 'uu', 'uu_xy', 'uu_xz', 'uu_yz' )
[3994]	969
[4039]	970	grid_x = 'xu'
	971	grid_y = 'y'
	972	grid_z = 'zu'
	973	!
	974	!-- v grid
	975	CASE ( 'vv', 'vv_xy', 'vv_xz', 'vv_yz' )
	976
	977	grid_x = 'x'
	978	grid_y = 'yv'
	979	grid_z = 'zu'
	980	!
	981	!-- w grid
[4583]	982	CASE ( 'ww', 'ww_xy', 'ww_xz', 'ww_yz', &
	983	'wu', 'wu_xy', 'wu_xz', 'wu_yz', &
	984	'wv', 'wv_xy', 'wv_xz', 'wv_yz', &
	985	'wtheta', 'wtheta_xy', 'wtheta_xz', 'wtheta_yz', &
[4866]	986	'wq', 'wq_xy', 'wq_xz', 'wq_yz', &
	987	'ws', 'ws_xy', 'ws_xz', 'ws_yz' )
[4039]	988
	989	grid_x = 'x'
	990	grid_y = 'y'
	991	grid_z = 'zw'
	992
	993	CASE DEFAULT
	994	found = .FALSE.
	995	grid_x = 'none'
	996	grid_y = 'none'
	997	grid_z = 'none'
	998
	999	END SELECT
	1000
	1001
	1002	END SUBROUTINE doq_define_netcdf_grid
[4431]	1003
[4583]	1004	!--------------------------------------------------------------------------------------------------!
[4039]	1005	!
	1006	! Description:
	1007	! ------------
	1008	!> Subroutine defining 2D output variables
[4583]	1009	!--------------------------------------------------------------------------------------------------!
	1010	SUBROUTINE doq_output_2d( av, variable, found, grid, mode, local_pf, two_d, nzb_do, nzt_do, &
	1011	fill_value )
[4039]	1012
	1013
	1014	IMPLICIT NONE
	1015
[4431]	1016	CHARACTER (LEN=*) :: grid !<
	1017	CHARACTER (LEN=*) :: mode !<
	1018	CHARACTER (LEN=*) :: variable !<
[4039]	1019
	1020	INTEGER(iwp) :: av !< value indicating averaged or non-averaged output
	1021	INTEGER(iwp) :: flag_nr !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
	1022	INTEGER(iwp) :: i !< grid index x-direction
	1023	INTEGER(iwp) :: j !< grid index y-direction
	1024	INTEGER(iwp) :: k !< grid index z-direction
[4431]	1025	INTEGER(iwp) :: nzb_do !<
	1026	INTEGER(iwp) :: nzt_do !<
[4039]	1027
	1028	LOGICAL :: found !< true if variable is in list
	1029	LOGICAL :: resorted !< true if array is resorted
	1030	LOGICAL :: two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
	1031
	1032	REAL(wp) :: fill_value !< value for the _FillValue attribute
[4431]	1033
	1034	REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !<
[4039]	1035	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which needs to be resorted for output
[4431]	1036
[4583]	1037
[4039]	1038	flag_nr = 0
	1039	found = .TRUE.
	1040	resorted = .FALSE.
	1041	two_d = .FALSE.
	1042
	1043	SELECT CASE ( TRIM( variable ) )
	1044
[4757]	1045	CASE ( 'rh_xy', 'rh_xz', 'rh_yz' )
	1046	IF ( av == 0 ) THEN
	1047	to_be_resorted => rh
	1048	ELSE
	1049	IF ( .NOT. ALLOCATED( rh_av ) ) THEN
	1050	ALLOCATE( rh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	1051	rh_av = REAL( fill_value, KIND = wp )
	1052	ENDIF
	1053	to_be_resorted => rh_av
	1054	ENDIF
	1055	flag_nr = 0
	1056
	1057	IF ( mode == 'xy' ) grid = 'zu'
	1058
[4861]	1059	CASE ( 'ta_xy', 'ta_xz', 'ta_yz' )
	1060	IF ( av == 0 ) THEN
	1061	to_be_resorted => ta
	1062	ELSE
	1063	IF ( .NOT. ALLOCATED( ta_av ) ) THEN
	1064	ALLOCATE( ta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	1065	ta_av = REAL( fill_value, KIND = wp )
	1066	ENDIF
	1067	to_be_resorted => ta_av
	1068	ENDIF
	1069	flag_nr = 0
	1070
	1071	IF ( mode == 'xy' ) grid = 'zu'
	1072
[4039]	1073	CASE ( 'ti_xy', 'ti_xz', 'ti_yz' )
	1074	IF ( av == 0 ) THEN
	1075	to_be_resorted => ti
	1076	ELSE
[4583]	1077	IF ( .NOT. ALLOCATED( ti_av ) ) THEN
[4518]	1078	ALLOCATE( ti_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4039]	1079	ti_av = REAL( fill_value, KIND = wp )
	1080	ENDIF
	1081	to_be_resorted => ti_av
	1082	ENDIF
	1083	flag_nr = 0
[4431]	1084
[4039]	1085	IF ( mode == 'xy' ) grid = 'zu'
[4431]	1086
[4039]	1087	CASE ( 'uu_xy', 'uu_xz', 'uu_yz' )
	1088	IF ( av == 0 ) THEN
	1089	to_be_resorted => uu
	1090	ELSE
[4583]	1091	IF ( .NOT. ALLOCATED( uu_av ) ) THEN
[4518]	1092	ALLOCATE( uu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4039]	1093	uu_av = REAL( fill_value, KIND = wp )
	1094	ENDIF
	1095	to_be_resorted => uu_av
	1096	ENDIF
	1097	flag_nr = 1
[4431]	1098
[4039]	1099	IF ( mode == 'xy' ) grid = 'zu'
[4431]	1100
[4039]	1101	CASE ( 'vv_xy', 'vv_xz', 'vv_yz' )
	1102	IF ( av == 0 ) THEN
	1103	to_be_resorted => vv
	1104	ELSE
[4583]	1105	IF ( .NOT. ALLOCATED( vv_av ) ) THEN
[4518]	1106	ALLOCATE( vv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4039]	1107	vv_av = REAL( fill_value, KIND = wp )
	1108	ENDIF
	1109	to_be_resorted => vv_av
	1110	ENDIF
	1111	flag_nr = 2
[4431]	1112
[4039]	1113	IF ( mode == 'xy' ) grid = 'zu'
[4431]	1114
[4039]	1115	CASE ( 'ww_xy', 'ww_xz', 'ww_yz' )
	1116	IF ( av == 0 ) THEN
	1117	to_be_resorted => ww
	1118	ELSE
[4583]	1119	IF ( .NOT. ALLOCATED( ww_av ) ) THEN
[4518]	1120	ALLOCATE( ww_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4039]	1121	ww_av = REAL( fill_value, KIND = wp )
	1122	ENDIF
	1123	to_be_resorted => ww_av
	1124	ENDIF
	1125	flag_nr = 3
[4431]	1126
[4039]	1127	IF ( mode == 'xy' ) grid = 'zw'
[4431]	1128
[4351]	1129	CASE ( 'wu_xy', 'wu_xz', 'wu_yz' )
	1130	IF ( av == 0 ) THEN
	1131	to_be_resorted => wu
	1132	ELSE
[4583]	1133	IF ( .NOT. ALLOCATED( wu_av ) ) THEN
[4518]	1134	ALLOCATE( wu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4351]	1135	wu_av = REAL( fill_value, KIND = wp )
	1136	ENDIF
	1137	to_be_resorted => wu_av
	1138	ENDIF
	1139	flag_nr = 0
[4431]	1140
[4351]	1141	IF ( mode == 'xy' ) grid = 'zw'
[4431]	1142
[4351]	1143	CASE ( 'wv_xy', 'wv_xz', 'wv_yz' )
	1144	IF ( av == 0 ) THEN
	1145	to_be_resorted => wv
	1146	ELSE
[4583]	1147	IF ( .NOT. ALLOCATED( wv_av ) ) THEN
[4518]	1148	ALLOCATE( wv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4351]	1149	wv_av = REAL( fill_value, KIND = wp )
	1150	ENDIF
	1151	to_be_resorted => wv_av
	1152	ENDIF
	1153	flag_nr = 0
[4431]	1154
[4351]	1155	IF ( mode == 'xy' ) grid = 'zw'
[4431]	1156
[4351]	1157	CASE ( 'wtheta_xy', 'wtheta_xz', 'wtheta_yz' )
	1158	IF ( av == 0 ) THEN
	1159	to_be_resorted => wtheta
	1160	ELSE
[4583]	1161	IF ( .NOT. ALLOCATED( wtheta_av ) ) THEN
[4518]	1162	ALLOCATE( wtheta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4351]	1163	wtheta_av = REAL( fill_value, KIND = wp )
	1164	ENDIF
	1165	to_be_resorted => wtheta_av
	1166	ENDIF
	1167	flag_nr = 0
[4431]	1168
[4351]	1169	IF ( mode == 'xy' ) grid = 'zw'
[4431]	1170
[4351]	1171	CASE ( 'wq_xy', 'wq_xz', 'wq_yz' )
	1172	IF ( av == 0 ) THEN
	1173	to_be_resorted => wq
	1174	ELSE
[4583]	1175	IF ( .NOT. ALLOCATED( wq_av ) ) THEN
[4518]	1176	ALLOCATE( wq_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4351]	1177	wq_av = REAL( fill_value, KIND = wp )
	1178	ENDIF
	1179	to_be_resorted => wq_av
	1180	ENDIF
	1181	flag_nr = 0
[4431]	1182
[4351]	1183	IF ( mode == 'xy' ) grid = 'zw'
[4431]	1184
[4866]	1185	CASE ( 'ws_xy', 'ws_xz', 'ws_yz' )
	1186	IF ( av == 0 ) THEN
	1187	to_be_resorted => ws
	1188	ELSE
	1189	IF ( .NOT. ALLOCATED( ws_av ) ) THEN
	1190	ALLOCATE( ws_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	1191	ws_av = REAL( fill_value, KIND = wp )
	1192	ENDIF
	1193	to_be_resorted => ws_av
	1194	ENDIF
	1195	flag_nr = 0
	1196
	1197	IF ( mode == 'xy' ) grid = 'zw'
	1198
[4331]	1199	CASE ( 'theta_2m*_xy' ) ! 2d-array
	1200	IF ( av == 0 ) THEN
	1201	DO i = nxl, nxr
	1202	DO j = nys, nyn
	1203	local_pf(i,j,nzb+1) = pt_2m(j,i)
	1204	ENDDO
	1205	ENDDO
	1206	ELSE
[4583]	1207	IF ( .NOT. ALLOCATED( pt_2m_av ) ) THEN
[4331]	1208	ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )
	1209	pt_2m_av = REAL( fill_value, KIND = wp )
	1210	ENDIF
	1211	DO i = nxl, nxr
	1212	DO j = nys, nyn
	1213	local_pf(i,j,nzb+1) = pt_2m_av(j,i)
	1214	ENDDO
	1215	ENDDO
	1216	ENDIF
	1217	resorted = .TRUE.
	1218	two_d = .TRUE.
	1219	grid = 'zu1'
[4431]	1220
[4331]	1221	CASE ( 'wspeed_10m*_xy' ) ! 2d-array
	1222	IF ( av == 0 ) THEN
	1223	DO i = nxl, nxr
	1224	DO j = nys, nyn
	1225	local_pf(i,j,nzb+1) = uv_10m(j,i)
	1226	ENDDO
	1227	ENDDO
	1228	ELSE
[4583]	1229	IF ( .NOT. ALLOCATED( uv_10m_av ) ) THEN
[4331]	1230	ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) )
	1231	uv_10m_av = REAL( fill_value, KIND = wp )
	1232	ENDIF
	1233	DO i = nxl, nxr
	1234	DO j = nys, nyn
	1235	local_pf(i,j,nzb+1) = uv_10m_av(j,i)
	1236	ENDDO
	1237	ENDDO
	1238	ENDIF
	1239	resorted = .TRUE.
	1240	two_d = .TRUE.
	1241	grid = 'zu1'
[4039]	1242
[4431]	1243	CASE ( 'wspeed_xy', 'wspeed_xz', 'wspeed_yz' )
	1244	IF ( av == 0 ) THEN
	1245	to_be_resorted => wspeed
	1246	ELSE
[4583]	1247	IF ( .NOT. ALLOCATED( wspeed_av ) ) THEN
[4518]	1248	ALLOCATE( wspeed_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4431]	1249	wspeed_av = REAL( fill_value, KIND = wp )
	1250	ENDIF
	1251	to_be_resorted => wspeed_av
	1252	ENDIF
	1253	flag_nr = 0
	1254
	1255	IF ( mode == 'xy' ) grid = 'zu'
	1256
	1257	CASE ( 'wdir_xy', 'wdir_xz', 'wdir_yz' )
	1258	IF ( av == 0 ) THEN
	1259	to_be_resorted => wdir
	1260	ELSE
[4583]	1261	IF ( .NOT. ALLOCATED( wdir_av ) ) THEN
[4518]	1262	ALLOCATE( wdir_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4431]	1263	wdir_av = REAL( fill_value, KIND = wp )
	1264	ENDIF
	1265	to_be_resorted => wdir_av
	1266	ENDIF
	1267	flag_nr = 0
	1268
	1269	IF ( mode == 'xy' ) grid = 'zu'
	1270
[4039]	1271	CASE DEFAULT
	1272	found = .FALSE.
	1273	grid = 'none'
	1274
	1275	END SELECT
[4431]	1276
	1277	IF ( found .AND. .NOT. resorted ) THEN
[4039]	1278	DO i = nxl, nxr
	1279	DO j = nys, nyn
	1280	DO k = nzb_do, nzt_do
[4583]	1281	local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), &
	1282	REAL( fill_value, KIND = wp ), &
	1283	BTEST( wall_flags_total_0(k,j,i), flag_nr ) )
[4039]	1284	ENDDO
	1285	ENDDO
	1286	ENDDO
[3994]	1287	ENDIF
[4431]	1288
[4039]	1289	END SUBROUTINE doq_output_2d
[4431]	1290
	1291
[4583]	1292	!--------------------------------------------------------------------------------------------------!
[4039]	1293	!
	1294	! Description:
	1295	! ------------
	1296	!> Subroutine defining 3D output variables
[4583]	1297	!--------------------------------------------------------------------------------------------------!
	1298	SUBROUTINE doq_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
[4431]	1299
[4039]	1300	IMPLICIT NONE
[3994]	1301
[4431]	1302	CHARACTER (LEN=*) :: variable !<
[3994]	1303
[4039]	1304	INTEGER(iwp) :: av !< index indicating averaged or instantaneous output
	1305	INTEGER(iwp) :: flag_nr !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
	1306	INTEGER(iwp) :: i !< index variable along x-direction
	1307	INTEGER(iwp) :: j !< index variable along y-direction
	1308	INTEGER(iwp) :: k !< index variable along z-direction
	1309	INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0)
	1310	INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d)
[3994]	1311
[4039]	1312	LOGICAL :: found !< true if variable is in list
	1313	LOGICAL :: resorted !< true if array is resorted
[3994]	1314
[4039]	1315	REAL(wp) :: fill_value !< value for the _FillValue attribute
	1316
[4768]	1317	REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !<
[4039]	1318	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which needs to be resorted for output
	1319
[4583]	1320
[4039]	1321	flag_nr = 0
	1322	found = .TRUE.
	1323	resorted = .FALSE.
[4431]	1324
[4039]	1325	SELECT CASE ( TRIM( variable ) )
	1326
[4757]	1327	CASE ( 'rh' )
	1328	IF ( av == 0 ) THEN
	1329	to_be_resorted => rh
	1330	ELSE
	1331	IF ( .NOT. ALLOCATED( rh ) ) THEN
	1332	ALLOCATE( rh(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	1333	rh = REAL( fill_value, KIND = wp )
	1334	ENDIF
	1335	to_be_resorted => rh
	1336	ENDIF
	1337	flag_nr = 0
	1338
[4861]	1339	CASE ( 'ta' )
	1340	IF ( av == 0 ) THEN
	1341	to_be_resorted => ta
	1342	ELSE
	1343	IF ( .NOT. ALLOCATED( ta ) ) THEN
	1344	ALLOCATE( ta(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	1345	ta = REAL( fill_value, KIND = wp )
	1346	ENDIF
	1347	to_be_resorted => ta
	1348	ENDIF
	1349	flag_nr = 0
	1350
[4039]	1351	CASE ( 'ti' )
	1352	IF ( av == 0 ) THEN
	1353	to_be_resorted => ti
	1354	ELSE
[4583]	1355	IF ( .NOT. ALLOCATED( ti_av ) ) THEN
[4518]	1356	ALLOCATE( ti_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4039]	1357	ti_av = REAL( fill_value, KIND = wp )
	1358	ENDIF
	1359	to_be_resorted => ti_av
	1360	ENDIF
	1361	flag_nr = 0
[4431]	1362
[4039]	1363	CASE ( 'uu' )
	1364	IF ( av == 0 ) THEN
	1365	to_be_resorted => uu
	1366	ELSE
[4583]	1367	IF ( .NOT. ALLOCATED( uu_av ) ) THEN
[4518]	1368	ALLOCATE( uu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4039]	1369	uu_av = REAL( fill_value, KIND = wp )
	1370	ENDIF
	1371	to_be_resorted => uu_av
	1372	ENDIF
	1373	flag_nr = 1
[4431]	1374
[4039]	1375	CASE ( 'vv' )
	1376	IF ( av == 0 ) THEN
	1377	to_be_resorted => vv
	1378	ELSE
[4583]	1379	IF ( .NOT. ALLOCATED( vv_av ) ) THEN
[4518]	1380	ALLOCATE( vv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4039]	1381	vv_av = REAL( fill_value, KIND = wp )
	1382	ENDIF
	1383	to_be_resorted => vv_av
	1384	ENDIF
	1385	flag_nr = 2
[4431]	1386
[4039]	1387	CASE ( 'ww' )
	1388	IF ( av == 0 ) THEN
	1389	to_be_resorted => ww
	1390	ELSE
[4583]	1391	IF ( .NOT. ALLOCATED( ww_av ) ) THEN
[4518]	1392	ALLOCATE( ww_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4039]	1393	ww_av = REAL( fill_value, KIND = wp )
	1394	ENDIF
	1395	to_be_resorted => ww_av
	1396	ENDIF
	1397	flag_nr = 3
	1398
[4351]	1399	CASE ( 'wu' )
	1400	IF ( av == 0 ) THEN
	1401	to_be_resorted => wu
	1402	ELSE
[4583]	1403	IF ( .NOT. ALLOCATED( wu_av ) ) THEN
[4518]	1404	ALLOCATE( wu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4351]	1405	wu_av = REAL( fill_value, KIND = wp )
	1406	ENDIF
	1407	to_be_resorted => wu_av
	1408	ENDIF
	1409	flag_nr = 0
	1410
	1411	CASE ( 'wv' )
	1412	IF ( av == 0 ) THEN
	1413	to_be_resorted => wv
	1414	ELSE
[4583]	1415	IF ( .NOT. ALLOCATED( wv_av ) ) THEN
[4518]	1416	ALLOCATE( wv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4351]	1417	wv_av = REAL( fill_value, KIND = wp )
	1418	ENDIF
	1419	to_be_resorted => wv_av
	1420	ENDIF
	1421	flag_nr = 0
[4431]	1422
[4351]	1423	CASE ( 'wtheta' )
	1424	IF ( av == 0 ) THEN
	1425	to_be_resorted => wtheta
	1426	ELSE
[4583]	1427	IF ( .NOT. ALLOCATED( wtheta_av ) ) THEN
[4518]	1428	ALLOCATE( wtheta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4351]	1429	wtheta_av = REAL( fill_value, KIND = wp )
	1430	ENDIF
	1431	to_be_resorted => wtheta_av
	1432	ENDIF
	1433	flag_nr = 0
[4431]	1434
[4351]	1435	CASE ( 'wq' )
	1436	IF ( av == 0 ) THEN
	1437	to_be_resorted => wq
	1438	ELSE
[4583]	1439	IF ( .NOT. ALLOCATED( wq_av ) ) THEN
[4518]	1440	ALLOCATE( wq_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4351]	1441	wq_av = REAL( fill_value, KIND = wp )
	1442	ENDIF
	1443	to_be_resorted => wq_av
	1444	ENDIF
	1445	flag_nr = 0
	1446
[4866]	1447	CASE ( 'ws' )
	1448	IF ( av == 0 ) THEN
	1449	to_be_resorted => ws
	1450	ELSE
	1451	IF ( .NOT. ALLOCATED( ws_av ) ) THEN
	1452	ALLOCATE( ws_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	1453	ws_av = REAL( fill_value, KIND = wp )
	1454	ENDIF
	1455	to_be_resorted => ws_av
	1456	ENDIF
	1457	flag_nr = 0
	1458
[4431]	1459	CASE ( 'wspeed' )
	1460	IF ( av == 0 ) THEN
	1461	to_be_resorted => wspeed
	1462	ELSE
[4583]	1463	IF ( .NOT. ALLOCATED( wspeed_av ) ) THEN
[4518]	1464	ALLOCATE( wspeed_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4431]	1465	wspeed_av = REAL( fill_value, KIND = wp )
	1466	ENDIF
	1467	to_be_resorted => wspeed_av
	1468	ENDIF
	1469	flag_nr = 0
	1470
	1471	CASE ( 'wdir' )
	1472	IF ( av == 0 ) THEN
	1473	to_be_resorted => wdir
	1474	ELSE
[4583]	1475	IF ( .NOT. ALLOCATED( wdir_av ) ) THEN
[4518]	1476	ALLOCATE( wdir_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4431]	1477	wdir_av = REAL( fill_value, KIND = wp )
	1478	ENDIF
	1479	to_be_resorted => wdir_av
	1480	ENDIF
	1481	flag_nr = 0
	1482
[4039]	1483	CASE DEFAULT
	1484	found = .FALSE.
	1485
	1486	END SELECT
[4431]	1487
	1488	IF ( found .AND. .NOT. resorted ) THEN
[4039]	1489	DO i = nxl, nxr
	1490	DO j = nys, nyn
	1491	DO k = nzb_do, nzt_do
[4583]	1492	local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), &
	1493	REAL( fill_value, KIND = wp ), &
	1494	BTEST( wall_flags_total_0(k,j,i), flag_nr ) )
[4039]	1495	ENDDO
	1496	ENDDO
	1497	ENDDO
	1498	ENDIF
	1499
	1500	END SUBROUTINE doq_output_3d
[4431]	1501
[4583]	1502
	1503	!--------------------------------------------------------------------------------------------------!
	1504	!
[3994]	1505	! Description:
	1506	! ------------
[4583]	1507	!> Resorts the user-defined output quantity with indices (k,j,i) to a temporary array with indices
	1508	!> (i,j,k) for masked data output.
	1509	!--------------------------------------------------------------------------------------------------!
[4069]	1510	SUBROUTINE doq_output_mask( av, variable, found, local_pf, mid )
[4431]	1511
[4039]	1512	USE control_parameters
[4431]	1513
[4039]	1514	USE indices
[4167]	1515
[4039]	1516	IMPLICIT NONE
[3994]	1517
[4583]	1518	REAL(wp), PARAMETER :: fill_value = -9999.0_wp !< value for the _FillValue attribute
	1519
[4167]	1520	CHARACTER (LEN=*) :: variable !<
	1521	CHARACTER (LEN=5) :: grid !< flag to distinquish between staggered grids
[3994]	1522
[4167]	1523	INTEGER(iwp) :: av !< index indicating averaged or instantaneous output
	1524	INTEGER(iwp) :: flag_nr !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
	1525	INTEGER(iwp) :: i !< index variable along x-direction
	1526	INTEGER(iwp) :: j !< index variable along y-direction
	1527	INTEGER(iwp) :: k !< index variable along z-direction
	1528	INTEGER(iwp) :: im !< loop index for masked variables
	1529	INTEGER(iwp) :: jm !< loop index for masked variables
	1530	INTEGER(iwp) :: kk !< masked output index variable along z-direction
	1531	INTEGER(iwp) :: mid !< masked output running index
	1532	INTEGER(iwp) :: ktt !< k index of highest horizontal surface
[3994]	1533
[4167]	1534	LOGICAL :: found !< true if variable is in list
	1535	LOGICAL :: resorted !< true if array is resorted
[3994]	1536
[4583]	1537	REAL(wp), DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) :: local_pf !<
[3994]	1538
[4583]	1539	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which needs to be resorted for output
[4167]	1540
[4583]	1541
[4039]	1542	flag_nr = 0
	1543	found = .TRUE.
	1544	resorted = .FALSE.
	1545	grid = 's'
	1546
	1547	SELECT CASE ( TRIM( variable ) )
	1548
[4757]	1549	CASE ( 'rh' )
	1550	IF ( av == 0 ) THEN
	1551	to_be_resorted => rh
	1552	ELSE
	1553	to_be_resorted => rh_av
	1554	ENDIF
	1555	grid = 's'
	1556	flag_nr = 0
	1557
[4861]	1558	CASE ( 'ta' )
	1559	IF ( av == 0 ) THEN
	1560	to_be_resorted => ta
	1561	ELSE
	1562	to_be_resorted => ta_av
	1563	ENDIF
	1564	grid = 's'
	1565	flag_nr = 0
	1566
[4039]	1567	CASE ( 'ti' )
	1568	IF ( av == 0 ) THEN
	1569	to_be_resorted => ti
	1570	ELSE
	1571	to_be_resorted => ti_av
	1572	ENDIF
	1573	grid = 's'
	1574	flag_nr = 0
[4431]	1575
[4039]	1576	CASE ( 'uu' )
	1577	IF ( av == 0 ) THEN
	1578	to_be_resorted => uu
	1579	ELSE
	1580	to_be_resorted => uu_av
	1581	ENDIF
	1582	grid = 'u'
	1583	flag_nr = 1
[4431]	1584
[4039]	1585	CASE ( 'vv' )
	1586	IF ( av == 0 ) THEN
	1587	to_be_resorted => vv
	1588	ELSE
	1589	to_be_resorted => vv_av
	1590	ENDIF
	1591	grid = 'v'
	1592	flag_nr = 2
[4431]	1593
[4039]	1594	CASE ( 'ww' )
	1595	IF ( av == 0 ) THEN
	1596	to_be_resorted => ww
	1597	ELSE
	1598	to_be_resorted => ww_av
	1599	ENDIF
	1600	grid = 'w'
	1601	flag_nr = 3
[4431]	1602
[4351]	1603	CASE ( 'wu' )
	1604	IF ( av == 0 ) THEN
	1605	to_be_resorted => wu
	1606	ELSE
	1607	to_be_resorted => wu_av
	1608	ENDIF
	1609	grid = 's'
	1610	flag_nr = 0
[4431]	1611
[4351]	1612	CASE ( 'wv' )
	1613	IF ( av == 0 ) THEN
	1614	to_be_resorted => wv
	1615	ELSE
	1616	to_be_resorted => wv_av
	1617	ENDIF
	1618	grid = 's'
	1619	flag_nr = 0
[4431]	1620
[4351]	1621	CASE ( 'wtheta' )
	1622	IF ( av == 0 ) THEN
	1623	to_be_resorted => wtheta
	1624	ELSE
	1625	to_be_resorted => wtheta_av
	1626	ENDIF
	1627	grid = 's'
	1628	flag_nr = 0
[4431]	1629
[4351]	1630	CASE ( 'wq' )
	1631	IF ( av == 0 ) THEN
	1632	to_be_resorted => wq
	1633	ELSE
	1634	to_be_resorted => wq_av
	1635	ENDIF
	1636	grid = 's'
	1637	flag_nr = 0
[4039]	1638
[4866]	1639	CASE ( 'ws' )
	1640	IF ( av == 0 ) THEN
	1641	to_be_resorted => ws
	1642	ELSE
	1643	to_be_resorted => ws_av
	1644	ENDIF
	1645	grid = 's'
	1646	flag_nr = 0
	1647
[4431]	1648	CASE ( 'wspeed' )
	1649	IF ( av == 0 ) THEN
	1650	to_be_resorted => wspeed
	1651	ELSE
	1652	to_be_resorted => wspeed_av
	1653	ENDIF
	1654	grid = 's'
	1655	flag_nr = 0
	1656
	1657	CASE ( 'wdir' )
	1658	IF ( av == 0 ) THEN
	1659	to_be_resorted => wdir
	1660	ELSE
	1661	to_be_resorted => wdir_av
	1662	ENDIF
	1663	grid = 's'
	1664	flag_nr = 0
	1665
[4039]	1666	CASE DEFAULT
	1667	found = .FALSE.
	1668
	1669	END SELECT
[4431]	1670
[4039]	1671	IF ( found .AND. .NOT. resorted ) THEN
	1672	IF ( .NOT. mask_surface(mid) ) THEN
	1673	!
	1674	!-- Default masked output
	1675	DO i = 1, mask_size_l(mid,1)
	1676	DO j = 1, mask_size_l(mid,2)
	1677	DO k = 1, mask_size_l(mid,3)
[4583]	1678	local_pf(i,j,k) = MERGE( to_be_resorted(mask_k(mid,k), &
	1679	mask_j(mid,j), &
	1680	mask_i(mid,i)), &
	1681	REAL( fill_value, KIND = wp ), &
	1682	BTEST( wall_flags_total_0(mask_k(mid,k), &
	1683	mask_j(mid,j), &
	1684	mask_i(mid,i)), &
[4431]	1685	flag_nr ) )
[4039]	1686	ENDDO
	1687	ENDDO
	1688	ENDDO
	1689
	1690	ELSE
	1691	!
	1692	!-- Terrain-following masked output
	1693	DO i = 1, mask_size_l(mid,1)
	1694	DO j = 1, mask_size_l(mid,2)
	1695	!
[4895]	1696	!-- Get k index of the lowest non-terrain grid point, else set fill value
[4167]	1697	im = mask_i(mid,i)
	1698	jm = mask_j(mid,j)
[4583]	1699	ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 ) ), DIM=1 ) - 1
[4167]	1700	DO k = 1, mask_size_l(mid,3)
[4895]	1701	kk = MIN( ktt + mask_k(mid,k) - 1, nzt+1 )
[4039]	1702	!
[4167]	1703	!-- Set value if not in building
[4346]	1704	IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) ) THEN
[4167]	1705	local_pf(i,j,k) = fill_value
	1706	ELSE
	1707	local_pf(i,j,k) = to_be_resorted(kk,jm,im)
	1708	ENDIF
[4039]	1709	ENDDO
	1710	ENDDO
	1711	ENDDO
	1712
	1713	ENDIF
	1714	ENDIF
[4431]	1715
[4039]	1716	END SUBROUTINE doq_output_mask
	1717
[4583]	1718	!--------------------------------------------------------------------------------------------------!
[4039]	1719	! Description:
	1720	! ------------
	1721	!> Allocate required arrays
[4583]	1722	!--------------------------------------------------------------------------------------------------!
[4039]	1723	SUBROUTINE doq_init
	1724
[3994]	1725	IMPLICIT NONE
[4431]	1726
[4039]	1727	INTEGER(iwp) :: ivar !< loop index over all 2d/3d/mask output quantities
[4157]	1728
[4583]	1729
[4039]	1730	!
	1731	!-- Next line is to avoid compiler warnings about unused variables
	1732	IF ( timestep_number_at_prev_calc == 0 ) CONTINUE
[4157]	1733	!
	1734	!-- Preparatory steps and initialization of output arrays
	1735	IF ( .NOT. prepared_diagnostic_output_quantities ) CALL doq_prepare
[3994]	1736
[4039]	1737	initialized_diagnostic_output_quantities = .FALSE.
[4431]	1738
[4039]	1739	ivar = 1
	1740
[4431]	1741	DO WHILE ( ivar <= SIZE( do_all ) )
	1742
[4039]	1743	SELECT CASE ( TRIM( do_all(ivar) ) )
	1744	!
[4757]	1745	!-- Allocate array for 'relative humidity'
	1746	CASE ( 'rh' )
	1747	IF ( .NOT. ALLOCATED( rh ) ) THEN
	1748	ALLOCATE( rh(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	1749	rh = 0.0_wp
	1750	ENDIF
	1751	!
[4861]	1752	!-- Allocate array for 'air temperature'
	1753	CASE ( 'ta' )
	1754	IF ( .NOT. ALLOCATED( ta ) ) THEN
	1755	ALLOCATE( ta(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	1756	ta = 0.0_wp
	1757	ENDIF
	1758	!
[4039]	1759	!-- Allocate array for 'turbulence intensity'
	1760	CASE ( 'ti' )
	1761	IF ( .NOT. ALLOCATED( ti ) ) THEN
[4518]	1762	ALLOCATE( ti(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4039]	1763	ti = 0.0_wp
	1764	ENDIF
	1765	!
	1766	!-- Allocate array for uu
	1767	CASE ( 'uu' )
	1768	IF ( .NOT. ALLOCATED( uu ) ) THEN
[4518]	1769	ALLOCATE( uu(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4039]	1770	uu = 0.0_wp
	1771	ENDIF
	1772	!
	1773	!-- Allocate array for vv
	1774	CASE ( 'vv' )
	1775	IF ( .NOT. ALLOCATED( vv ) ) THEN
[4518]	1776	ALLOCATE( vv(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4039]	1777	vv = 0.0_wp
	1778	ENDIF
	1779	!
	1780	!-- Allocate array for ww
	1781	CASE ( 'ww' )
	1782	IF ( .NOT. ALLOCATED( ww ) ) THEN
[4518]	1783	ALLOCATE( ww(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4039]	1784	ww = 0.0_wp
	1785	ENDIF
[4331]	1786	!
[4351]	1787	!-- Allocate array for wu
	1788	CASE ( 'wu' )
	1789	IF ( .NOT. ALLOCATED( wu ) ) THEN
[4518]	1790	ALLOCATE( wu(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4351]	1791	wu = 0.0_wp
	1792	ENDIF
	1793	!
	1794	!-- Allocate array for wv
	1795	CASE ( 'wv' )
	1796	IF ( .NOT. ALLOCATED( wv ) ) THEN
[4518]	1797	ALLOCATE( wv(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4351]	1798	wv = 0.0_wp
	1799	ENDIF
	1800	!
	1801	!-- Allocate array for wtheta
	1802	CASE ( 'wtheta' )
	1803	IF ( .NOT. ALLOCATED( wtheta ) ) THEN
[4518]	1804	ALLOCATE( wtheta(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4351]	1805	wtheta = 0.0_wp
	1806	ENDIF
	1807	!
	1808	!-- Allocate array for wq
	1809	CASE ( 'wq' )
	1810	IF ( .NOT. ALLOCATED( wq ) ) THEN
[4518]	1811	ALLOCATE( wq(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4351]	1812	wq = 0.0_wp
	1813	ENDIF
	1814	!
[4866]	1815	!-- Allocate array for ws
	1816	CASE ( 'ws' )
	1817	IF ( .NOT. ALLOCATED( ws ) ) THEN
	1818	ALLOCATE( ws(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	1819	ws = 0.0_wp
	1820	ENDIF
	1821	!
[4331]	1822	!-- Allocate array for 2-m potential temperature
	1823	CASE ( 'theta_2m*' )
	1824	IF ( .NOT. ALLOCATED( pt_2m ) ) THEN
[4518]	1825	ALLOCATE( pt_2m(nysg:nyng,nxlg:nxrg) )
[4331]	1826	pt_2m = 0.0_wp
	1827	ENDIF
	1828	!
	1829	!-- Allocate array for 10-m wind speed
	1830	CASE ( 'wspeed_10m*' )
	1831	IF ( .NOT. ALLOCATED( uv_10m ) ) THEN
[4518]	1832	ALLOCATE( uv_10m(nysg:nyng,nxlg:nxrg) )
[4331]	1833	uv_10m = 0.0_wp
	1834	ENDIF
[4431]	1835	!
	1836	!-- Allocate array for wspeed
	1837	CASE ( 'wspeed' )
	1838	IF ( .NOT. ALLOCATED( wspeed ) ) THEN
[4518]	1839	ALLOCATE( wspeed(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4431]	1840	wspeed = 0.0_wp
	1841	ENDIF
[4039]	1842
[4431]	1843	!
	1844	!-- Allocate array for wdir
	1845	CASE ( 'wdir' )
	1846	IF ( .NOT. ALLOCATED( u_center ) ) THEN
[4518]	1847	ALLOCATE( u_center(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4431]	1848	u_center = 0.0_wp
	1849	ENDIF
	1850	IF ( .NOT. ALLOCATED( v_center ) ) THEN
[4518]	1851	ALLOCATE( v_center(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4431]	1852	v_center = 0.0_wp
	1853	ENDIF
	1854	IF ( .NOT. ALLOCATED( wdir ) ) THEN
[4518]	1855	ALLOCATE( wdir(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[4431]	1856	wdir = 0.0_wp
	1857	ENDIF
	1858
[4039]	1859	END SELECT
	1860
	1861	ivar = ivar + 1
	1862	ENDDO
	1863
	1864	initialized_diagnostic_output_quantities = .TRUE.
[4431]	1865
[4039]	1866	END SUBROUTINE doq_init
	1867
	1868
	1869	!--------------------------------------------------------------------------------------------------!
	1870	! Description:
	1871	! ------------
[4757]	1872	!> Calculation of module-specific statistics, i.e. horizontally averaged profiles and time series.
	1873	!> This is called for every statistic region sr, but at least for the region "total domain" (sr=0).
	1874	!--------------------------------------------------------------------------------------------------!
	1875	SUBROUTINE doq_statistics( mode, sr, tn )
	1876
	1877
	1878	CHARACTER (LEN=*) :: mode !<
	1879
	1880	! INTEGER(iwp) :: i !<
	1881	! INTEGER(iwp) :: j !<
	1882	! INTEGER(iwp) :: k !<
	1883	INTEGER(iwp) :: sr !<
	1884	INTEGER(iwp) :: tn !<
	1885	!
	1886	!-- Next line is to avoid compiler warning about unused variables. Please remove.
	1887	IF ( sr == 0 .OR. tn == 0 ) CONTINUE
	1888
	1889	IF ( mode == 'profiles' ) THEN
	1890
	1891	ELSEIF ( mode == 'time_series' ) THEN
	1892
	1893	ENDIF
	1894
	1895	END SUBROUTINE doq_statistics
	1896
	1897
	1898	!--------------------------------------------------------------------------------------------------!
	1899	! Description:
	1900	! ------------
[4039]	1901	!> Calculate diagnostic quantities
	1902	!--------------------------------------------------------------------------------------------------!
	1903	SUBROUTINE doq_calculate
	1904
	1905	IMPLICIT NONE
	1906
[4331]	1907	INTEGER(iwp) :: i !< grid index x-dimension
[4431]	1908	INTEGER(iwp) :: j !< grid index y-dimension
[4331]	1909	INTEGER(iwp) :: k !< grid index z-dimension
[4039]	1910	INTEGER(iwp) :: ivar !< loop index over all 2d/3d/mask output quantities
[4431]	1911
[4757]	1912	REAL(wp) :: temp !< temperature
	1913	REAL(wp) :: e_s !< saturation vapor pressure
	1914	REAL(wp) :: q_s !< saturation mixing ratio
	1915
[4331]	1916	TYPE(surf_type), POINTER :: surf !< surf-type array, used to generalize subroutines
[3994]	1917
	1918
	1919	! CALL cpu_log( log_point(41), 'calculate_quantities', 'start' )
[4157]	1920
[3994]	1921	!
[4583]	1922	!-- Save timestep number to check in time_integration if doq_calculate has been called already,
	1923	!-- since the CALL occurs at two locations, but the calculations need to be done only once per
	1924	!-- timestep.
[3994]	1925	timestep_number_at_prev_calc = current_timestep_number
	1926
[4039]	1927	ivar = 1
[3994]	1928
[4431]	1929	DO WHILE ( ivar <= SIZE( do_all ) )
[3994]	1930
[4039]	1931	SELECT CASE ( TRIM( do_all(ivar) ) )
[3994]	1932	!
[4757]	1933	!-- rh (relative humidity)
	1934	CASE ( 'rh' )
	1935	DO i = nxl, nxr
	1936	DO j = nys, nyn
	1937	DO k = nzb+1, nzt
	1938	IF ( humidity ) THEN
	1939	temp = exner(k) * pt(k,j,i)
	1940	e_s = magnus( temp )
	1941	q_s = rd_d_rv * e_s / ( hyp(k) - e_s )
	1942	rh(k,j,i) = q(k,j,i) / q_s * 100.0_wp &
	1943	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
	1944	ENDIF
	1945	ENDDO
	1946	ENDDO
	1947	ENDDO
	1948	!
[4861]	1949	!-- ta (air temperature)
	1950	CASE ( 'ta' )
	1951	DO i = nxl, nxr
	1952	DO j = nys, nyn
	1953	DO k = nzb+1, nzt
	1954	ta(k,j,i) = exner(k) * pt(k,j,i) - degc_to_k
	1955	ENDDO
	1956	ENDDO
	1957	ENDDO
	1958	!
[3994]	1959	!-- Calculate 'turbulence intensity' from rot[(u,v,w)] at scalar grid point
	1960	CASE ( 'ti' )
	1961	DO i = nxl, nxr
	1962	DO j = nys, nyn
	1963	DO k = nzb+1, nzt
[4583]	1964	ti(k,j,i) = 0.25_wp * SQRT( &
	1965	( ( w(k,j+1,i) + w(k-1,j+1,i) &
	1966	- w(k,j-1,i) - w(k-1,j-1,i) ) * ddy &
	1967	- ( v(k+1,j,i) + v(k+1,j+1,i) &
	1968	- v(k-1,j,i) - v(k-1,j+1,i) ) * ddzu(k) )**2 &
	1969	+ ( ( u(k+1,j,i) + u(k+1,j,i+1) &
	1970	- u(k-1,j,i) - u(k-1,j,i+1) ) * ddzu(k) &
	1971	- ( w(k,j,i+1) + w(k-1,j,i+1) &
	1972	- w(k,j,i-1) - w(k-1,j,i-1) ) * ddx )**2 &
	1973	+ ( ( v(k,j,i+1) + v(k,j+1,i+1) &
	1974	- v(k,j,i-1) - v(k,j+1,i-1) ) * ddx &
	1975	- ( u(k,j+1,i) + u(k,j+1,i+1) &
	1976	- u(k,j-1,i) - u(k,j-1,i+1) ) * ddy )**2 ) &
	1977	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
[3994]	1978	ENDDO
	1979	ENDDO
[4431]	1980	ENDDO
[4039]	1981	!
	1982	!-- uu
	1983	CASE ( 'uu' )
	1984	DO i = nxl, nxr
	1985	DO j = nys, nyn
	1986	DO k = nzb+1, nzt
[4583]	1987	uu(k,j,i) = u(k,j,i) * u(k,j,i) &
	1988	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
[4039]	1989	ENDDO
	1990	ENDDO
[3994]	1991	ENDDO
[4039]	1992	!
	1993	!-- vv
	1994	CASE ( 'vv' )
	1995	DO i = nxl, nxr
	1996	DO j = nys, nyn
	1997	DO k = nzb+1, nzt
[4583]	1998	vv(k,j,i) = v(k,j,i) * v(k,j,i) &
	1999	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
[4039]	2000	ENDDO
	2001	ENDDO
	2002	ENDDO
	2003	!
	2004	!-- ww
	2005	CASE ( 'ww' )
	2006	DO i = nxl, nxr
	2007	DO j = nys, nyn
	2008	DO k = nzb+1, nzt-1
[4583]	2009	ww(k,j,i) = w(k,j,i) * w(k,j,i) &
	2010	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
[4039]	2011	ENDDO
	2012	ENDDO
	2013	ENDDO
[4331]	2014	!
[4351]	2015	!-- wu
	2016	CASE ( 'wu' )
	2017	DO i = nxl, nxr
	2018	DO j = nys, nyn
	2019	DO k = nzb+1, nzt-1
[4583]	2020	wu(k,j,i) = w(k,j,i) &
	2021	* 0.25_wp * ( u(k,j,i) + u(k,j,i+1) + u(k+1,j,i) + u(k+1,j,i+1) )&
	2022	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
[4351]	2023	ENDDO
	2024	ENDDO
	2025	ENDDO
	2026	!
	2027	!-- wv
	2028	CASE ( 'wv' )
	2029	DO i = nxl, nxr
	2030	DO j = nys, nyn
	2031	DO k = nzb+1, nzt-1
[4583]	2032	wv(k,j,i) = w(k,j,i) &
	2033	* 0.25_wp * ( v(k,j,i) + v(k,j+1,i) + v(k+1,j,i) + v(k+1,j+1,i) )&
	2034	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
[4351]	2035	ENDDO
	2036	ENDDO
	2037	ENDDO
	2038	!
	2039	!-- wtheta
	2040	CASE ( 'wtheta' )
	2041	DO i = nxl, nxr
	2042	DO j = nys, nyn
	2043	DO k = nzb+1, nzt-1
[4583]	2044	wtheta(k,j,i) = w(k,j,i) &
	2045	* 0.5_wp * ( pt(k,j,i) + pt(k+1,j,i) ) &
	2046	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ))
[4351]	2047	ENDDO
	2048	ENDDO
	2049	ENDDO
	2050	!
	2051	!-- wq
	2052	CASE ( 'wq' )
	2053	DO i = nxl, nxr
	2054	DO j = nys, nyn
	2055	DO k = nzb+1, nzt-1
[4583]	2056	wq(k,j,i) = w(k,j,i) * 0.5_wp * ( q(k,j,i) + q(k+1,j,i) ) &
	2057	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
[4351]	2058	ENDDO
	2059	ENDDO
	2060	ENDDO
	2061	!
[4866]	2062	!-- ws
	2063	CASE ( 'ws' )
	2064	DO i = nxl, nxr
	2065	DO j = nys, nyn
	2066	DO k = nzb+1, nzt-1
	2067	ws(k,j,i) = w(k,j,i) * 0.5_wp * ( s(k,j,i) + s(k+1,j,i) ) &
	2068	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
	2069	ENDDO
	2070	ENDDO
	2071	ENDDO
	2072	!
[4331]	2073	!-- 2-m potential temperature
	2074	CASE ( 'theta_2m*' )
	2075	!
[4583]	2076	!-- 2-m potential temperature is caluclated from surface arrays. In case the 2m level is
	2077	!-- below the first grid point, MOST is applied, else, linear interpolation between two
	2078	!-- vertical grid levels is applied. To access all surfaces, iterate over all horizontally-
[4331]	2079	!-- upward facing surface types.
	2080	surf => surf_def_h(0)
	2081	CALL calc_pt_2m
[4671]	2082	surf => surf_lsm_h(0)
[4331]	2083	CALL calc_pt_2m
[4671]	2084	surf => surf_usm_h(0)
[4331]	2085	CALL calc_pt_2m
	2086	!
	2087	!-- 10-m wind speed
	2088	CASE ( 'wspeed_10m*' )
	2089	!
[4583]	2090	!-- 10-m wind speed is caluclated from surface arrays. In case the 10m level is below the
	2091	!-- first grid point, MOST is applied, else, linear interpolation between two vertical grid
	2092	!-- levels is applied. To access all surfaces, iterate over all horizontally-upward facing
	2093	!-- surface types.
[4331]	2094	surf => surf_def_h(0)
	2095	CALL calc_wind_10m
[4671]	2096	surf => surf_lsm_h(0)
[4331]	2097	CALL calc_wind_10m
[4671]	2098	surf => surf_usm_h(0)
[4331]	2099	CALL calc_wind_10m
[4431]	2100	!
[4583]	2101	!-- Horizontal wind speed
[4431]	2102	CASE ( 'wspeed' )
	2103	DO i = nxl, nxr
	2104	DO j = nys, nyn
	2105	DO k = nzb, nzt+1
[4583]	2106	wspeed(k,j,i) = SQRT( ( 0.5_wp * ( u(k,j,i) + u(k,j,i+1) ) )**2 &
	2107	+ ( 0.5_wp * ( v(k,j,i) + v(k,j+1,i) ) )**2 ) &
	2108	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ))
[4431]	2109	ENDDO
	2110	ENDDO
	2111	ENDDO
[3994]	2112
[4431]	2113	!
[4583]	2114	!-- Horizontal wind direction
[4431]	2115	CASE ( 'wdir' )
	2116	DO i = nxl, nxr
	2117	DO j = nys, nyn
	2118	DO k = nzb, nzt+1
	2119	u_center(k,j,i) = 0.5_wp * ( u(k,j,i) + u(k,j,i+1) )
	2120	v_center(k,j,i) = 0.5_wp * ( v(k,j,i) + v(k,j+1,i) )
	2121
[4583]	2122	wdir(k,j,i) = ATAN2( u_center(k,j,i), v_center(k,j,i) ) &
	2123	/ pi * 180.0_wp + 180.0_wp
[4431]	2124	ENDDO
	2125	ENDDO
	2126	ENDDO
	2127
[4039]	2128	END SELECT
[3994]	2129
[4039]	2130	ivar = ivar + 1
[3994]	2131	ENDDO
	2132
	2133	! CALL cpu_log( log_point(41), 'calculate_quantities', 'stop' )
	2134
[4331]	2135	!
[4583]	2136	!-- The following block contains subroutines to calculate diagnostic surface quantities.
[4331]	2137	CONTAINS
[4583]	2138	!--------------------------------------------------------------------------------------------------!
[4331]	2139	! Description:
	2140	! ------------
	2141	!> Calculation of 2-m potential temperature.
[4583]	2142	!--------------------------------------------------------------------------------------------------!
[4331]	2143	SUBROUTINE calc_pt_2m
	2144
[4583]	2145	USE surface_layer_fluxes_mod, &
[4331]	2146	ONLY: psi_h
	2147
	2148	IMPLICIT NONE
	2149
	2150	INTEGER(iwp) :: kk !< running index along the z-dimension
	2151	INTEGER(iwp) :: m !< running index for surface elements
[4431]	2152
[4331]	2153	DO m = 1, surf%ns
	2154
	2155	i = surf%i(m)
	2156	j = surf%j(m)
	2157	k = surf%k(m)
	2158	!
[4583]	2159	!-- If 2-m level is below the first grid level, MOST is used for calculation of
	2160	!-- 2-m temperature.
[4331]	2161	IF ( surf%z_mo(m) > 2.0_wp ) THEN
[4583]	2162	pt_2m(j,i) = surf%pt_surface(m) + surf%ts(m) / kappa &
	2163	* ( LOG( 2.0_wp / surf%z0h(m) ) &
	2164	- psi_h( 2.0_wp / surf%ol(m) ) &
	2165	+ psi_h( surf%z0h(m) / surf%ol(m) ) )
[4331]	2166	!
[4583]	2167	!-- If 2-m level is above the first grid level, 2-m temperature is linearly interpolated
	2168	!-- between the two nearest vertical grid levels. Note, since 2-m temperature is only
	2169	!-- computed for horizontal upward-facing surfaces, only a vertical interpolation is
	2170	!-- necessary.
[4331]	2171	ELSE
	2172	!
[4431]	2173	!-- zw(k-1) defines the height of the surface.
[4331]	2174	kk = k
	2175	DO WHILE ( zu(kk) - zw(k-1) < 2.0_wp .AND. kk <= nzt )
[4431]	2176	kk = kk + 1
[4331]	2177	ENDDO
	2178	!
[4431]	2179	!-- kk defines the index of the first grid level >= 2m.
[4583]	2180	pt_2m(j,i) = pt(kk-1,j,i) + &
	2181	( zw(k-1) + 2.0_wp - zu(kk-1) ) * &
	2182	( pt(kk,j,i) - pt(kk-1,j,i) ) / &
[4331]	2183	( zu(kk) - zu(kk-1) )
	2184	ENDIF
	2185
	2186	ENDDO
	2187
	2188	END SUBROUTINE calc_pt_2m
[4431]	2189
[4583]	2190	!--------------------------------------------------------------------------------------------------!
[4331]	2191	! Description:
	2192	! ------------
	2193	!> Calculation of 10-m wind speed.
[4583]	2194	!--------------------------------------------------------------------------------------------------!
[4331]	2195	SUBROUTINE calc_wind_10m
	2196
[4583]	2197	USE surface_layer_fluxes_mod, &
[4331]	2198	ONLY: psi_m
	2199
	2200	IMPLICIT NONE
	2201
	2202	INTEGER(iwp) :: kk !< running index along the z-dimension
	2203	INTEGER(iwp) :: m !< running index for surface elements
	2204
	2205	REAL(wp) :: uv_l !< wind speed at lower grid point
	2206	REAL(wp) :: uv_u !< wind speed at upper grid point
[4431]	2207
[4583]	2208
[4331]	2209	DO m = 1, surf%ns
	2210
	2211	i = surf%i(m)
	2212	j = surf%j(m)
	2213	k = surf%k(m)
	2214	!
[4583]	2215	!-- If 10-m level is below the first grid level, MOST is used for calculation of 10-m
	2216	!-- temperature.
[4331]	2217	IF ( surf%z_mo(m) > 10.0_wp ) THEN
[4583]	2218	uv_10m(j,i) = surf%us(m) / kappa &
	2219	* ( LOG( 10.0_wp / surf%z0(m) ) &
	2220	- psi_m( 10.0_wp / surf%ol(m) ) &
	2221	+ psi_m( surf%z0(m) / surf%ol(m) ) )
[4331]	2222	!
[4583]	2223	!-- If 10-m level is above the first grid level, 10-m wind speed is linearly interpolated
	2224	!-- between the two nearest vertical grid levels. Note, since 10-m temperature is only
	2225	!-- computed for horizontal upward-facing surfaces, only a vertical interpolation is
	2226	!-- necessary.
[4331]	2227	ELSE
	2228	!
[4431]	2229	!-- zw(k-1) defines the height of the surface.
[4331]	2230	kk = k
	2231	DO WHILE ( zu(kk) - zw(k-1) < 10.0_wp .AND. kk <= nzt )
[4431]	2232	kk = kk + 1
[4331]	2233	ENDDO
	2234	!
	2235	!-- kk defines the index of the first grid level >= 10m.
[4583]	2236	uv_l = SQRT( ( 0.5_wp * ( u(kk-1,j,i) + u(kk-1,j,i+1) ) )**2 &
[4331]	2237	+ ( 0.5_wp * ( v(kk-1,j,i) + v(kk-1,j+1,i) ) )**2 )
	2238
[4583]	2239	uv_u = SQRT( ( 0.5_wp * ( u(kk,j,i) + u(kk,j,i+1) ) )**2 &
[4331]	2240	+ ( 0.5_wp * ( v(kk,j,i) + v(kk,j+1,i) ) )**2 )
	2241
[4583]	2242	uv_10m(j,i) = uv_l + ( zw(k-1) + 10.0_wp - zu(kk-1) ) * &
	2243	( uv_u - uv_l ) / &
[4331]	2244	( zu(kk) - zu(kk-1) )
	2245
	2246	ENDIF
	2247
	2248	ENDDO
	2249
	2250	END SUBROUTINE calc_wind_10m
	2251
[4039]	2252	END SUBROUTINE doq_calculate
[3994]	2253
	2254
[4583]	2255	!--------------------------------------------------------------------------------------------------!
[3994]	2256	! Description:
	2257	! ------------
[4583]	2258	!> Preparation of the diagnostic output, counting of the module-specific output quantities and
	2259	!> gathering of the output names.
	2260	!--------------------------------------------------------------------------------------------------!
[4039]	2261	SUBROUTINE doq_prepare
[3994]	2262
[4583]	2263	USE control_parameters, &
[4069]	2264	ONLY: do2d, do3d, domask, masks
[3994]	2265
	2266	IMPLICIT NONE
	2267
[4583]	2268	CHARACTER (LEN=varnamelength), DIMENSION(0:1,500) :: do2d_var = ' ' !< label array for 2d output quantities
[3994]	2269
	2270	INTEGER(iwp) :: av !< index defining type of output, av=0 instantaneous, av=1 averaged
	2271	INTEGER(iwp) :: ivar !< loop index
	2272	INTEGER(iwp) :: ivar_all !< loop index
	2273	INTEGER(iwp) :: l !< index for cutting string
[4431]	2274	INTEGER(iwp) :: mid !< masked output running index
[3994]	2275
[4583]	2276
[3994]	2277	prepared_diagnostic_output_quantities = .FALSE.
	2278
	2279	ivar = 1
	2280	ivar_all = 1
	2281
	2282	DO av = 0, 1
	2283	!
	2284	!-- Remove _xy, _xz, or _yz from string
	2285	l = MAX( 3, LEN_TRIM( do2d(av,ivar) ) )
[3998]	2286	do2d_var(av,ivar)(1:l-3) = do2d(av,ivar)(1:l-3)
[3994]	2287	!
[4039]	2288	!-- Gather 2d output quantity names.
[4583]	2289	!-- Check for double occurrence of output quantity, e.g. by _xy, _yz, _xz.
[3994]	2290	DO WHILE ( do2d_var(av,ivar)(1:1) /= ' ' )
	2291	IF ( .NOT. ANY( do_all == do2d_var(av,ivar) ) ) THEN
	2292	do_all(ivar_all) = do2d_var(av,ivar)
	2293	ENDIF
	2294	ivar = ivar + 1
	2295	ivar_all = ivar_all + 1
	2296	l = MAX( 3, LEN_TRIM( do2d(av,ivar) ) )
[3998]	2297	do2d_var(av,ivar)(1:l-3) = do2d(av,ivar)(1:l-3)
[3994]	2298	ENDDO
	2299
	2300	ivar = 1
	2301	!
[4039]	2302	!-- Gather 3d output quantity names
[3994]	2303	DO WHILE ( do3d(av,ivar)(1:1) /= ' ' )
	2304	do_all(ivar_all) = do3d(av,ivar)
	2305	ivar = ivar + 1
	2306	ivar_all = ivar_all + 1
	2307	ENDDO
	2308
	2309	ivar = 1
	2310	!
[4583]	2311	!-- Gather masked output quantity names. Also check for double output e.g. by different masks.
[3994]	2312	DO mid = 1, masks
	2313	DO WHILE ( domask(mid,av,ivar)(1:1) /= ' ' )
[4039]	2314	IF ( .NOT. ANY( do_all == domask(mid,av,ivar) ) ) THEN
[4132]	2315	do_all(ivar_all) = domask(mid,av,ivar)
[4039]	2316	ENDIF
[3994]	2317
	2318	ivar = ivar + 1
	2319	ivar_all = ivar_all + 1
	2320	ENDDO
	2321	ivar = 1
	2322	ENDDO
	2323
	2324	ENDDO
	2325
	2326	prepared_diagnostic_output_quantities = .TRUE.
	2327
[4039]	2328	END SUBROUTINE doq_prepare
[4431]	2329
[4583]	2330	!--------------------------------------------------------------------------------------------------!
[4039]	2331	! Description:
	2332	! ------------
	2333	!> Subroutine reads local (subdomain) restart data
[4583]	2334	!--------------------------------------------------------------------------------------------------!
	2335	SUBROUTINE doq_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, nxr_on_file, nynf, nync, &
	2336	nyn_on_file, nysf, nysc, nys_on_file, tmp_2d, tmp_3d, found )
[4431]	2337
[4518]	2338	USE control_parameters
	2339
	2340	USE indices
	2341
	2342	USE kinds
	2343
	2344	USE pegrid
	2345
	2346
	2347	IMPLICIT NONE
	2348
	2349	INTEGER(iwp) :: k !<
	2350	INTEGER(iwp) :: nxlc !<
	2351	INTEGER(iwp) :: nxlf !<
	2352	INTEGER(iwp) :: nxl_on_file !<
	2353	INTEGER(iwp) :: nxrc !<
	2354	INTEGER(iwp) :: nxrf !<
	2355	INTEGER(iwp) :: nxr_on_file !<
	2356	INTEGER(iwp) :: nync !<
	2357	INTEGER(iwp) :: nynf !<
	2358	INTEGER(iwp) :: nyn_on_file !<
	2359	INTEGER(iwp) :: nysc !<
	2360	INTEGER(iwp) :: nysf !<
	2361	INTEGER(iwp) :: nys_on_file !<
	2362
	2363	LOGICAL, INTENT(OUT) :: found
	2364
	2365	REAL(wp), DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_2d !<
	2366
	2367	REAL(wp), DIMENSION(nzb:nzt+1,nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: tmp_3d !<
	2368
	2369
	2370	found = .TRUE.
	2371
	2372	SELECT CASE ( restart_string(1:length) )
	2373
	2374	CASE ( 'pt_2m_av' )
	2375	IF ( .NOT. ALLOCATED( pt_2m_av ) ) THEN
	2376	ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )
	2377	ENDIF
	2378	IF ( k == 1 ) READ ( 13 ) tmp_2d
	2379	pt_2m_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	2380	tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
[4583]	2381
[4757]	2382	CASE ( 'rh_av' )
	2383	IF ( .NOT. ALLOCATED( rh_av ) ) THEN
	2384	ALLOCATE( rh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2385	ENDIF
	2386	IF ( k == 1 ) READ ( 13 ) tmp_3d
	2387	rh_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	2388	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	2389
[4861]	2390	CASE ( 'ta_av' )
	2391	IF ( .NOT. ALLOCATED( ta_av ) ) THEN
	2392	ALLOCATE( ta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2393	ENDIF
	2394	IF ( k == 1 ) READ ( 13 ) tmp_3d
	2395	ta_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	2396	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	2397
[4518]	2398	CASE ( 'ti_av' )
	2399	IF ( .NOT. ALLOCATED( ti_av ) ) THEN
	2400	ALLOCATE( ti_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2401	ENDIF
	2402	IF ( k == 1 ) READ ( 13 ) tmp_3d
	2403	ti_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	2404	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	2405
	2406	CASE ( 'u_center_av' )
	2407	IF ( .NOT. ALLOCATED( u_center_av ) ) THEN
	2408	ALLOCATE( u_center_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2409	ENDIF
	2410	IF ( k == 1 ) READ ( 13 ) tmp_3d
	2411	u_center_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	2412	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	2413
	2414	CASE ( 'uu_av' )
	2415	IF ( .NOT. ALLOCATED( uu_av ) ) THEN
	2416	ALLOCATE( uu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2417	ENDIF
	2418	IF ( k == 1 ) READ ( 13 ) tmp_3d
	2419	uu_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	2420	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	2421
	2422	CASE ( 'uv_10m_av' )
	2423	IF ( .NOT. ALLOCATED( uv_10m_av ) ) THEN
	2424	ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) )
	2425	ENDIF
	2426	IF ( k == 1 ) READ ( 13 ) tmp_2d
	2427	uv_10m_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	2428	tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	2429
	2430	CASE ( 'v_center_av' )
	2431	IF ( .NOT. ALLOCATED( v_center_av ) ) THEN
	2432	ALLOCATE( v_center_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2433	ENDIF
	2434	IF ( k == 1 ) READ ( 13 ) tmp_3d
	2435	v_center_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	2436	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	2437
	2438	CASE ( 'vv_av' )
	2439	IF ( .NOT. ALLOCATED( vv_av ) ) THEN
	2440	ALLOCATE( vv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2441	ENDIF
	2442	IF ( k == 1 ) READ ( 13 ) tmp_3d
	2443	vv_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	2444	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	2445
	2446	CASE ( 'wtheta_av' )
	2447	IF ( .NOT. ALLOCATED( wtheta_av ) ) THEN
	2448	ALLOCATE( wtheta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2449	ENDIF
	2450	IF ( k == 1 ) READ ( 13 ) tmp_3d
	2451	wtheta_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	2452	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	2453
	2454	CASE ( 'wq_av' )
	2455	IF ( .NOT. ALLOCATED( wq_av ) ) THEN
	2456	ALLOCATE( wq_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2457	ENDIF
	2458	IF ( k == 1 ) READ ( 13 ) tmp_3d
	2459	wq_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	2460	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	2461
[4866]	2462	CASE ( 'ws_av' )
	2463	IF ( .NOT. ALLOCATED( ws_av ) ) THEN
	2464	ALLOCATE( ws_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2465	ENDIF
	2466	IF ( k == 1 ) READ ( 13 ) tmp_3d
	2467	ws_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	2468	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	2469
[4518]	2470	CASE ( 'wspeed_av' )
	2471	IF ( .NOT. ALLOCATED( wspeed_av ) ) THEN
	2472	ALLOCATE( wspeed_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2473	ENDIF
	2474	IF ( k == 1 ) READ ( 13 ) tmp_3d
	2475	wspeed_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	2476	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	2477
	2478	CASE ( 'wu_av' )
	2479	IF ( .NOT. ALLOCATED( wu_av ) ) THEN
	2480	ALLOCATE( wu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2481	ENDIF
	2482	IF ( k == 1 ) READ ( 13 ) tmp_3d
	2483	wu_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	2484	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	2485
	2486	CASE ( 'wv_av' )
	2487	IF ( .NOT. ALLOCATED( wv_av ) ) THEN
	2488	ALLOCATE( wv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2489	ENDIF
	2490	IF ( k == 1 ) READ ( 13 ) tmp_3d
	2491	wv_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	2492	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	2493
	2494	CASE ( 'ww_av' )
	2495	IF ( .NOT. ALLOCATED( ww_av ) ) THEN
	2496	ALLOCATE( ww_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2497	ENDIF
	2498	IF ( k == 1 ) READ ( 13 ) tmp_3d
	2499	ww_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
	2500	tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
	2501
	2502
	2503	CASE DEFAULT
	2504
	2505	found = .FALSE.
	2506
	2507	END SELECT
	2508
	2509	END SUBROUTINE doq_rrd_local_ftn
	2510
[4583]	2511	!--------------------------------------------------------------------------------------------------!
[4039]	2512	! Description:
	2513	! ------------
[4518]	2514	!> Read module-specific local restart data arrays (MPI-IO).
[4583]	2515	!--------------------------------------------------------------------------------------------------!
[4518]	2516	SUBROUTINE doq_rrd_local_mpi
	2517
	2518	LOGICAL :: array_found !< control flad indicating whether the array is found in the file or not
	2519
	2520
	2521	CALL rd_mpi_io_check_array( 'pt_2m_av' , found = array_found )
	2522	IF ( array_found ) THEN
	2523	IF ( .NOT. ALLOCATED( pt_2m_av ) ) ALLOCATE( pt_2m_av(nysg:nyng,nxlg:nxrg) )
	2524	CALL rrd_mpi_io( 'pt_2m_av', pt_2m_av )
	2525	ENDIF
	2526
[4861]	2527	CALL rd_mpi_io_check_array( 'ta_av' , found = array_found )
	2528	IF ( array_found ) THEN
	2529	IF ( .NOT. ALLOCATED( ta_av ) ) ALLOCATE( ta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2530	CALL rrd_mpi_io( 'ta_av', ta_av )
	2531	ENDIF
	2532
[4518]	2533	CALL rd_mpi_io_check_array( 'ti_av' , found = array_found )
	2534	IF ( array_found ) THEN
	2535	IF ( .NOT. ALLOCATED( ti_av ) ) ALLOCATE( ti_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2536	CALL rrd_mpi_io( 'ti_av', ti_av )
	2537	ENDIF
	2538
[4757]	2539	CALL rd_mpi_io_check_array( 'rh_av' , found = array_found )
	2540	IF ( array_found ) THEN
	2541	IF ( .NOT. ALLOCATED( rh_av ) ) ALLOCATE( ti_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2542	CALL rrd_mpi_io( 'rh_av', rh_av )
	2543	ENDIF
	2544
[4518]	2545	CALL rd_mpi_io_check_array( 'u_center_av' , found = array_found )
	2546	IF ( array_found ) THEN
	2547	IF ( .NOT. ALLOCATED( u_center_av ) ) ALLOCATE( u_center_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2548	CALL rrd_mpi_io( 'u_center_av', u_center_av )
	2549	ENDIF
	2550
	2551	CALL rd_mpi_io_check_array( 'uu_av' , found = array_found )
	2552	IF ( array_found ) THEN
	2553	IF ( .NOT. ALLOCATED( uu_av ) ) ALLOCATE( uu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2554	CALL rrd_mpi_io( 'uu_av', uu_av )
	2555	ENDIF
	2556
	2557	CALL rd_mpi_io_check_array( 'uv_10m_av' , found = array_found )
	2558	IF ( array_found ) THEN
	2559	IF ( .NOT. ALLOCATED( uv_10m_av ) ) ALLOCATE( uv_10m_av(nysg:nyng,nxlg:nxrg) )
	2560	CALL rrd_mpi_io( 'uv_10m_av', uv_10m_av )
	2561	ENDIF
	2562
	2563	CALL rd_mpi_io_check_array( 'v_center_av' , found = array_found )
	2564	IF ( array_found ) THEN
	2565	IF ( .NOT. ALLOCATED( v_center_av ) ) ALLOCATE( v_center_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2566	CALL rrd_mpi_io( 'v_center_av', v_center_av )
	2567	ENDIF
	2568
	2569	CALL rd_mpi_io_check_array( 'vv_av' , found = array_found )
	2570	IF ( array_found ) THEN
	2571	IF ( .NOT. ALLOCATED( vv_av ) ) ALLOCATE( vv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2572	CALL rrd_mpi_io( 'vv_av', vv_av )
	2573	ENDIF
	2574
	2575	CALL rd_mpi_io_check_array( 'ww_av' , found = array_found )
	2576	IF ( array_found ) THEN
	2577	IF ( .NOT. ALLOCATED( ww_av ) ) ALLOCATE( ww_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2578	CALL rrd_mpi_io( 'ww_av', ww_av )
	2579	ENDIF
	2580
	2581	CALL rd_mpi_io_check_array( 'wu_av' , found = array_found )
	2582	IF ( array_found ) THEN
	2583	IF ( .NOT. ALLOCATED( wu_av ) ) ALLOCATE( wu_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2584	CALL rrd_mpi_io( 'wu_av', wu_av )
	2585	ENDIF
	2586
	2587	CALL rd_mpi_io_check_array( 'wv_av' , found = array_found )
	2588	IF ( array_found ) THEN
	2589	IF ( .NOT. ALLOCATED( wv_av ) ) ALLOCATE( wv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2590	CALL rrd_mpi_io( 'wv_av', wv_av )
	2591	ENDIF
	2592
	2593	CALL rd_mpi_io_check_array( 'wtheta_av' , found = array_found )
	2594	IF ( array_found ) THEN
	2595	IF ( .NOT. ALLOCATED( wtheta_av ) ) ALLOCATE( wtheta_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2596	CALL rrd_mpi_io( 'wtheta_av', wtheta_av )
	2597	ENDIF
	2598
	2599	CALL rd_mpi_io_check_array( 'wq_av' , found = array_found )
	2600	IF ( array_found ) THEN
	2601	IF ( .NOT. ALLOCATED( wq_av ) ) ALLOCATE( wq_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2602	CALL rrd_mpi_io( 'wq_av', wq_av )
	2603	ENDIF
	2604
[4866]	2605	CALL rd_mpi_io_check_array( 'ws_av' , found = array_found )
	2606	IF ( array_found ) THEN
	2607	IF ( .NOT. ALLOCATED( ws_av ) ) ALLOCATE( ws_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2608	CALL rrd_mpi_io( 'ws_av', ws_av )
	2609	ENDIF
	2610
[4518]	2611	CALL rd_mpi_io_check_array( 'wspeed_av' , found = array_found )
	2612	IF ( array_found ) THEN
	2613	IF ( .NOT. ALLOCATED( wspeed_av ) ) ALLOCATE( wspeed_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
	2614	CALL rrd_mpi_io( 'wspeed_av', wspeed_av )
	2615	ENDIF
	2616
	2617	END SUBROUTINE doq_rrd_local_mpi
	2618
[4583]	2619	!--------------------------------------------------------------------------------------------------!
[4518]	2620	! Description:
	2621	! ------------
[4039]	2622	!> Subroutine writes local (subdomain) restart data
[4583]	2623	!--------------------------------------------------------------------------------------------------!
[4039]	2624	SUBROUTINE doq_wrd_local
[3994]	2625
[4039]	2626	IMPLICIT NONE
	2627
[4583]	2628
[4517]	2629	IF ( TRIM( restart_data_format_output ) == 'fortran_binary' ) THEN
[4331]	2630
[4517]	2631	IF ( ALLOCATED( pt_2m_av ) ) THEN
	2632	CALL wrd_write_string( 'pt_2m_av' )
	2633	WRITE ( 14 ) pt_2m_av
	2634	ENDIF
[4331]	2635
[4757]	2636	IF ( ALLOCATED( rh_av ) ) THEN
	2637	CALL wrd_write_string( 'rh_av' )
	2638	WRITE ( 14 ) rh_av
	2639	ENDIF
	2640
[4861]	2641	IF ( ALLOCATED( ta_av ) ) THEN
	2642	CALL wrd_write_string( 'ta_av' )
	2643	WRITE ( 14 ) ta_av
	2644	ENDIF
	2645
[4517]	2646	IF ( ALLOCATED( ti_av ) ) THEN
	2647	CALL wrd_write_string( 'ti_av' )
	2648	WRITE ( 14 ) ti_av
	2649	ENDIF
[4331]	2650
[4518]	2651	IF ( ALLOCATED( u_center_av ) ) THEN
	2652	CALL wrd_write_string( 'u_center_av' )
	2653	WRITE ( 14 ) u_center_av
	2654	ENDIF
	2655
[4517]	2656	IF ( ALLOCATED( uu_av ) ) THEN
	2657	CALL wrd_write_string( 'uu_av' )
	2658	WRITE ( 14 ) uu_av
	2659	ENDIF
[4331]	2660
[4517]	2661	IF ( ALLOCATED( uv_10m_av ) ) THEN
	2662	CALL wrd_write_string( 'uv_10m_av' )
	2663	WRITE ( 14 ) uv_10m_av
	2664	ENDIF
[4331]	2665
[4518]	2666	IF ( ALLOCATED( v_center_av ) ) THEN
	2667	CALL wrd_write_string( 'v_center_av' )
	2668	WRITE ( 14 ) v_center_av
	2669	ENDIF
	2670
[4517]	2671	IF ( ALLOCATED( vv_av ) ) THEN
	2672	CALL wrd_write_string( 'vv_av' )
	2673	WRITE ( 14 ) vv_av
	2674	ENDIF
[4039]	2675
[4517]	2676	IF ( ALLOCATED( ww_av ) ) THEN
	2677	CALL wrd_write_string( 'ww_av' )
	2678	WRITE ( 14 ) ww_av
	2679	ENDIF
[4039]	2680
[4517]	2681	IF ( ALLOCATED( wu_av ) ) THEN
	2682	CALL wrd_write_string( 'wu_av' )
	2683	WRITE ( 14 ) wu_av
	2684	ENDIF
[4431]	2685
[4517]	2686	IF ( ALLOCATED( wv_av ) ) THEN
	2687	CALL wrd_write_string( 'wv_av' )
	2688	WRITE ( 14 ) wv_av
	2689	ENDIF
[4431]	2690
[4517]	2691	IF ( ALLOCATED( wtheta_av ) ) THEN
	2692	CALL wrd_write_string( 'wtheta_av' )
	2693	WRITE ( 14 ) wtheta_av
	2694	ENDIF
[4351]	2695
[4517]	2696	IF ( ALLOCATED( wq_av ) ) THEN
	2697	CALL wrd_write_string( 'wq_av' )
	2698	WRITE ( 14 ) wq_av
	2699	ENDIF
[4431]	2700
[4866]	2701	IF ( ALLOCATED( ws_av ) ) THEN
	2702	CALL wrd_write_string( 'ws_av' )
	2703	WRITE ( 14 ) ws_av
	2704	ENDIF
	2705
[4517]	2706	IF ( ALLOCATED( wspeed_av ) ) THEN
	2707	CALL wrd_write_string( 'wspeed_av' )
	2708	WRITE ( 14 ) wspeed_av
	2709	ENDIF
[4431]	2710
[4535]	2711	ELSEIF ( restart_data_format_output(1:3) == 'mpi' ) THEN
[4517]	2712
[4518]	2713	IF ( ALLOCATED( pt_2m_av ) ) CALL wrd_mpi_io( 'pt_2m_av', pt_2m_av )
[4757]	2714	IF ( ALLOCATED( rh_av ) ) CALL wrd_mpi_io( 'rh_av', rh_av )
[4861]	2715	IF ( ALLOCATED( ta_av ) ) CALL wrd_mpi_io( 'ta_av', ta_av )
[4518]	2716	IF ( ALLOCATED( ti_av ) ) CALL wrd_mpi_io( 'ti_av', ti_av )
	2717	IF ( ALLOCATED( u_center_av ) ) CALL wrd_mpi_io( 'u_center_av', u_center_av )
	2718	IF ( ALLOCATED( uu_av ) ) CALL wrd_mpi_io( 'uu_av', uu_av )
	2719	IF ( ALLOCATED( uv_10m_av ) ) CALL wrd_mpi_io( 'uv_10m_av', uv_10m_av )
	2720	IF ( ALLOCATED( vv_av ) ) CALL wrd_mpi_io( 'vv_av', vv_av )
	2721	IF ( ALLOCATED( v_center_av ) ) CALL wrd_mpi_io( 'v_center_av', v_center_av )
	2722	IF ( ALLOCATED( wtheta_av ) ) CALL wrd_mpi_io( 'wtheta_av', wtheta_av )
	2723	IF ( ALLOCATED( wq_av ) ) CALL wrd_mpi_io( 'wq_av', wq_av )
[4866]	2724	IF ( ALLOCATED( ws_av ) ) CALL wrd_mpi_io( 'ws_av', ws_av )
[4518]	2725	IF ( ALLOCATED( wspeed_av ) ) CALL wrd_mpi_io( 'wspeed_av', wspeed_av )
	2726	IF ( ALLOCATED( wu_av ) ) CALL wrd_mpi_io( 'wu_av', wu_av )
	2727	IF ( ALLOCATED( wv_av ) ) CALL wrd_mpi_io( 'wv_av', wv_av )
	2728	IF ( ALLOCATED( ww_av ) ) CALL wrd_mpi_io( 'ww_av', ww_av )
[4517]	2729
[4431]	2730	ENDIF
	2731
[4039]	2732	END SUBROUTINE doq_wrd_local
	2733
[4431]	2734
[3994]	2735	END MODULE diagnostic_output_quantities_mod

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

Context Navigation

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

Download in other formats: