Home

Context Navigation

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

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