Home

Context Navigation

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

Last change on this file since 4427 was 4360, checked in by suehring, 5 years ago
Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update
Property svn:keywords set to `Id`
File size: 63.2 KB

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