Home

Context Navigation

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

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