Home

Context Navigation

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

Last change on this file since 4753 was 4671, checked in by pavelkrc, 4 years ago
Radiative transfer model RTM version 4.1
Property svn:keywords set to `Id`
File size: 81.2 KB

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