Home

Context Navigation

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

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