Home

Context Navigation

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

Last change on this file since 4517 was 4517, checked in by raasch, 4 years ago
added restart with MPI-IO for reading local arrays
Property svn:keywords set to `Id`
File size: 71.5 KB

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