Home

Context Navigation

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

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