Home

Context Navigation

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

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