Home

Context Navigation

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

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