Home

Context Navigation

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

Last change on this file since 4125 was 4069, checked in by Giersch, 5 years ago
Bugfix for masked output, compiler warning removed, test case for wind turbine model revised
Property svn:keywords set to `Id`
File size: 38.1 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 4069 2019-07-01 14:05:51Z suehring $
[4069]	27	! Masked output running index mid has been introduced as a local variable to
	28	! avoid runtime error (Loop variable has been modified) in time_integration
	29	!
	30	! 4039 2019-06-18 10:32:41Z suehring
[4039]	31	! - Add output of uu, vv, ww to enable variance calculation according temporal
	32	! EC method
	33	! - Allocate arrays only when they are required
	34	! - Formatting adjustment
	35	! - Rename subroutines
	36	! - Further modularization
	37	!
	38	! 3998 2019-05-23 13:38:11Z suehring
[3998]	39	! Bugfix in gathering all output strings
	40	!
	41	! 3995 2019-05-22 18:59:54Z suehring
[3995]	42	! Avoid compiler warnings about unused variable and fix string operation which
	43	! is not allowed with PGI compiler
	44	!
	45	! 3994 2019-05-22 18:08:09Z suehring
[3994]	46	! Initial revision
	47	!
	48	!
	49	! @author Farah Kanani-Suehring
	50	!
	51	! Description:
	52	! ------------
	53	!> ...
	54	!------------------------------------------------------------------------------!
	55	MODULE diagnostic_output_quantities_mod
	56
[4039]	57	USE arrays_3d, &
	58	ONLY: ddzu, u, v, w
[3994]	59	!
	60	! USE averaging
	61	!
	62	! USE basic_constants_and_equations_mod, &
	63	! ONLY: c_p, lv_d_cp, l_v
	64	!
	65	! USE bulk_cloud_model_mod, &
	66	! ONLY: bulk_cloud_model
	67	!
	68	USE control_parameters, &
	69	ONLY: current_timestep_number, varnamelength
	70	!
	71	! USE cpulog, &
	72	! ONLY: cpu_log, log_point
[4039]	73
	74	USE grid_variables, &
	75	ONLY: ddx, ddy
[3994]	76	!
[4039]	77	USE indices, &
	78	ONLY: nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
[3994]	79	!
	80	USE kinds
	81	!
	82	! USE land_surface_model_mod, &
	83	! ONLY: zs
	84	!
	85	! USE module_interface, &
	86	! ONLY: module_interface_data_output_2d
	87	!
	88	! #if defined( __netcdf )
	89	! USE NETCDF
	90	! #endif
	91	!
	92	! USE netcdf_interface, &
	93	! ONLY: fill_value, id_set_xy, id_set_xz, id_set_yz, id_var_do2d, &
	94	! id_var_time_xy, id_var_time_xz, id_var_time_yz, nc_stat, &
	95	! netcdf_data_format, netcdf_handle_error
	96	!
	97	! USE particle_attributes, &
	98	! ONLY: grid_particles, number_of_particles, particle_advection_start, &
	99	! particles, prt_count
	100	!
	101	! USE pegrid
	102	!
	103	! USE surface_mod, &
	104	! ONLY: ind_pav_green, ind_veg_wall, ind_wat_win, surf_def_h, &
	105	! surf_lsm_h, surf_usm_h
	106	!
	107	! USE turbulence_closure_mod, &
	108	! ONLY: tcm_data_output_2d
	109
	110
	111	IMPLICIT NONE
	112
	113	CHARACTER(LEN=varnamelength), DIMENSION(500) :: do_all = ' '
	114
	115	INTEGER(iwp) :: timestep_number_at_prev_calc = 0 !< ...at previous diagnostic output calculation
[4039]	116
	117	LOGICAL :: initialized_diagnostic_output_quantities = .FALSE. !< flag indicating whether output is initialized
	118	LOGICAL :: prepared_diagnostic_output_quantities = .FALSE. !< flag indicating whether output is p
[3994]	119
	120	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ti !< rotation(u,v,w) aka turbulence intensity
	121	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ti_av !< avg. rotation(u,v,w) aka turbulence intensity
[4039]	122
	123	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: uu !< uu
	124	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: uu_av !< mean of uu
	125
	126	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: vv !< vv
	127	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: vv_av !< mean of vv
	128
	129	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ww !< ww
	130	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ww_av !< mean of ww
[3994]	131
	132
	133	SAVE
	134
	135	PRIVATE
	136
	137	!
	138	!-- Public variables
[4039]	139	PUBLIC do_all, &
	140	initialized_diagnostic_output_quantities, &
	141	prepared_diagnostic_output_quantities, &
	142	timestep_number_at_prev_calc, &
	143	ti_av, &
	144	uu_av, &
	145	vv_av, &
	146	ww_av
	147	!
	148	!-- Public routines
	149	PUBLIC doq_3d_data_averaging, &
	150	doq_calculate, &
	151	doq_check_data_output, &
	152	doq_define_netcdf_grid, &
	153	doq_output_2d, &
	154	doq_output_3d, &
	155	doq_output_mask, &
	156	doq_init, &
	157	doq_prepare, &
	158	doq_wrd_local
	159	! doq_rrd_local, &
[3994]	160
	161
[4039]	162	INTERFACE doq_3d_data_averaging
	163	MODULE PROCEDURE doq_3d_data_averaging
	164	END INTERFACE doq_3d_data_averaging
[3994]	165
[4039]	166	INTERFACE doq_calculate
	167	MODULE PROCEDURE doq_calculate
	168	END INTERFACE doq_calculate
[3994]	169
[4039]	170	INTERFACE doq_check_data_output
	171	MODULE PROCEDURE doq_check_data_output
	172	END INTERFACE doq_check_data_output
	173
	174	INTERFACE doq_define_netcdf_grid
	175	MODULE PROCEDURE doq_define_netcdf_grid
	176	END INTERFACE doq_define_netcdf_grid
	177
	178	INTERFACE doq_output_2d
	179	MODULE PROCEDURE doq_output_2d
	180	END INTERFACE doq_output_2d
	181
	182	INTERFACE doq_output_3d
	183	MODULE PROCEDURE doq_output_3d
	184	END INTERFACE doq_output_3d
	185
	186	INTERFACE doq_output_mask
	187	MODULE PROCEDURE doq_output_mask
	188	END INTERFACE doq_output_mask
	189
	190	INTERFACE doq_init
	191	MODULE PROCEDURE doq_init
	192	END INTERFACE doq_init
[3994]	193
[4039]	194	INTERFACE doq_prepare
	195	MODULE PROCEDURE doq_prepare
	196	END INTERFACE doq_prepare
	197
	198	! INTERFACE doq_rrd_local
	199	! MODULE PROCEDURE doq_rrd_local
	200	! END INTERFACE doq_rrd_local
	201
	202	INTERFACE doq_wrd_local
	203	MODULE PROCEDURE doq_wrd_local
	204	END INTERFACE doq_wrd_local
[3994]	205
[4039]	206
[3994]	207	CONTAINS
[4039]	208
	209	!------------------------------------------------------------------------------!
	210	! Description:
	211	! ------------
	212	!> Sum up and time-average diagnostic output quantities as well as allocate
	213	!> the array necessary for storing the average.
	214	!------------------------------------------------------------------------------!
	215	SUBROUTINE doq_3d_data_averaging( mode, variable )
[3994]	216
[4039]	217	USE control_parameters, &
	218	ONLY: average_count_3d
	219
	220	CHARACTER (LEN=*) :: mode !<
	221	CHARACTER (LEN=*) :: variable !<
	222
	223	INTEGER(iwp) :: i !<
	224	INTEGER(iwp) :: j !<
	225	INTEGER(iwp) :: k !<
	226
	227	IF ( mode == 'allocate' ) THEN
	228
	229	SELECT CASE ( TRIM( variable ) )
	230
	231	CASE ( 'ti' )
	232	IF ( .NOT. ALLOCATED( ti_av ) ) THEN
	233	ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
	234	ENDIF
	235	ti_av = 0.0_wp
	236
	237	CASE ( 'uu' )
	238	IF ( .NOT. ALLOCATED( uu_av ) ) THEN
	239	ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
	240	ENDIF
	241	uu_av = 0.0_wp
	242
	243	CASE ( 'vv' )
	244	IF ( .NOT. ALLOCATED( vv_av ) ) THEN
	245	ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
	246	ENDIF
	247	vv_av = 0.0_wp
	248
	249	CASE ( 'ww' )
	250	IF ( .NOT. ALLOCATED( ww_av ) ) THEN
	251	ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
	252	ENDIF
	253	ww_av = 0.0_wp
	254
	255	CASE DEFAULT
	256	CONTINUE
	257
	258	END SELECT
	259
	260	ELSEIF ( mode == 'sum' ) THEN
	261
	262	SELECT CASE ( TRIM( variable ) )
	263
	264	CASE ( 'ti' )
	265	IF ( ALLOCATED( ti_av ) ) THEN
	266	DO i = nxl, nxr
	267	DO j = nys, nyn
	268	DO k = nzb, nzt+1
	269	ti_av(k,j,i) = ti_av(k,j,i) + ti(k,j,i)
	270	ENDDO
	271	ENDDO
	272	ENDDO
	273	ENDIF
	274
	275	CASE ( 'uu' )
	276	IF ( ALLOCATED( uu_av ) ) THEN
	277	DO i = nxl, nxr
	278	DO j = nys, nyn
	279	DO k = nzb, nzt+1
	280	uu_av(k,j,i) = uu_av(k,j,i) + uu(k,j,i)
	281	ENDDO
	282	ENDDO
	283	ENDDO
	284	ENDIF
	285
	286	CASE ( 'vv' )
	287	IF ( ALLOCATED( vv_av ) ) THEN
	288	DO i = nxl, nxr
	289	DO j = nys, nyn
	290	DO k = nzb, nzt+1
	291	vv_av(k,j,i) = vv_av(k,j,i) + vv(k,j,i)
	292	ENDDO
	293	ENDDO
	294	ENDDO
	295	ENDIF
	296
	297	CASE ( 'ww' )
	298	IF ( ALLOCATED( ww_av ) ) THEN
	299	DO i = nxl, nxr
	300	DO j = nys, nyn
	301	DO k = nzb, nzt+1
	302	ww_av(k,j,i) = ww_av(k,j,i) + ww(k,j,i)
	303	ENDDO
	304	ENDDO
	305	ENDDO
	306	ENDIF
	307
	308	CASE DEFAULT
	309	CONTINUE
	310
	311	END SELECT
	312
	313	ELSEIF ( mode == 'average' ) THEN
	314
	315	SELECT CASE ( TRIM( variable ) )
	316
	317	CASE ( 'ti' )
	318	IF ( ALLOCATED( ti_av ) ) THEN
	319	DO i = nxl, nxr
	320	DO j = nys, nyn
	321	DO k = nzb, nzt+1
	322	ti_av(k,j,i) = ti_av(k,j,i) / REAL( average_count_3d, KIND=wp )
	323	ENDDO
	324	ENDDO
	325	ENDDO
	326	ENDIF
	327
	328	CASE ( 'uu' )
	329	IF ( ALLOCATED( uu_av ) ) THEN
	330	DO i = nxl, nxr
	331	DO j = nys, nyn
	332	DO k = nzb, nzt+1
	333	uu_av(k,j,i) = uu_av(k,j,i) / REAL( average_count_3d, KIND=wp )
	334	ENDDO
	335	ENDDO
	336	ENDDO
	337	ENDIF
	338
	339	CASE ( 'vv' )
	340	IF ( ALLOCATED( vv_av ) ) THEN
	341	DO i = nxl, nxr
	342	DO j = nys, nyn
	343	DO k = nzb, nzt+1
	344	vv_av(k,j,i) = vv_av(k,j,i) / REAL( average_count_3d, KIND=wp )
	345	ENDDO
	346	ENDDO
	347	ENDDO
	348	ENDIF
	349
	350	CASE ( 'ww' )
	351	IF ( ALLOCATED( ww_av ) ) THEN
	352	DO i = nxl, nxr
	353	DO j = nys, nyn
	354	DO k = nzb, nzt+1
	355	ww_av(k,j,i) = ww_av(k,j,i) / REAL( average_count_3d, KIND=wp )
	356	ENDDO
	357	ENDDO
	358	ENDDO
	359	ENDIF
	360
	361	END SELECT
	362
	363	ENDIF
	364
	365
	366	END SUBROUTINE doq_3d_data_averaging
	367
[3994]	368	!------------------------------------------------------------------------------!
	369	! Description:
	370	! ------------
[4039]	371	!> Check data output for diagnostic output
[3994]	372	!------------------------------------------------------------------------------!
[4039]	373	SUBROUTINE doq_check_data_output( var, unit )
[3994]	374
[4039]	375	IMPLICIT NONE
[3994]	376
[4039]	377	CHARACTER (LEN=*) :: unit !<
	378	CHARACTER (LEN=*) :: var !<
[3994]	379
[4039]	380	SELECT CASE ( TRIM( var ) )
	381
	382	CASE ( 'ti' )
	383	unit = '1/s'
	384
	385	CASE ( 'uu' )
	386	unit = 'm2/s2'
	387
	388	CASE ( 'vv' )
	389	unit = 'm2/s2'
	390
	391	CASE ( 'ww' )
	392	unit = 'm2/s2'
	393
	394	CASE DEFAULT
	395	unit = 'illegal'
	396
	397	END SELECT
	398
	399
	400	END SUBROUTINE doq_check_data_output
	401
	402	!------------------------------------------------------------------------------!
	403	!
	404	! Description:
	405	! ------------
	406	!> Subroutine defining appropriate grid for netcdf variables.
	407	!------------------------------------------------------------------------------!
	408	SUBROUTINE doq_define_netcdf_grid( variable, found, grid_x, grid_y, grid_z )
	409
[3994]	410	IMPLICIT NONE
[4039]	411
	412	CHARACTER (LEN=*), INTENT(IN) :: variable !<
	413	LOGICAL, INTENT(OUT) :: found !<
	414	CHARACTER (LEN=*), INTENT(OUT) :: grid_x !<
	415	CHARACTER (LEN=*), INTENT(OUT) :: grid_y !<
	416	CHARACTER (LEN=*), INTENT(OUT) :: grid_z !<
	417
	418	found = .TRUE.
	419
	420	SELECT CASE ( TRIM( variable ) )
[3995]	421	!
[4039]	422	!-- s grid
	423	CASE ( 'ti', 'ti_xy', 'ti_xz', 'ti_yz' )
[3994]	424
[4039]	425	grid_x = 'x'
	426	grid_y = 'y'
	427	grid_z = 'zu'
	428	!
	429	!-- u grid
	430	CASE ( 'uu', 'uu_xy', 'uu_xz', 'uu_yz' )
[3994]	431
[4039]	432	grid_x = 'xu'
	433	grid_y = 'y'
	434	grid_z = 'zu'
	435	!
	436	!-- v grid
	437	CASE ( 'vv', 'vv_xy', 'vv_xz', 'vv_yz' )
	438
	439	grid_x = 'x'
	440	grid_y = 'yv'
	441	grid_z = 'zu'
	442	!
	443	!-- w grid
	444	CASE ( 'ww', 'ww_xy', 'ww_xz', 'ww_yz' )
	445
	446	grid_x = 'x'
	447	grid_y = 'y'
	448	grid_z = 'zw'
	449
	450	CASE DEFAULT
	451	found = .FALSE.
	452	grid_x = 'none'
	453	grid_y = 'none'
	454	grid_z = 'none'
	455
	456	END SELECT
	457
	458
	459	END SUBROUTINE doq_define_netcdf_grid
	460
	461	!------------------------------------------------------------------------------!
	462	!
	463	! Description:
	464	! ------------
	465	!> Subroutine defining 2D output variables
	466	!------------------------------------------------------------------------------!
	467	SUBROUTINE doq_output_2d( av, variable, found, grid, &
	468	mode, local_pf, two_d, nzb_do, nzt_do, fill_value )
	469
	470
	471	IMPLICIT NONE
	472
	473	CHARACTER (LEN=*) :: grid !<
	474	CHARACTER (LEN=*) :: mode !<
	475	CHARACTER (LEN=*) :: variable !<
	476
	477	INTEGER(iwp) :: av !< value indicating averaged or non-averaged output
	478	INTEGER(iwp) :: flag_nr !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
	479	INTEGER(iwp) :: i !< grid index x-direction
	480	INTEGER(iwp) :: j !< grid index y-direction
	481	INTEGER(iwp) :: k !< grid index z-direction
	482	INTEGER(iwp) :: nzb_do !<
	483	INTEGER(iwp) :: nzt_do !<
	484
	485	LOGICAL :: found !< true if variable is in list
	486	LOGICAL :: resorted !< true if array is resorted
	487	LOGICAL :: two_d !< flag parameter that indicates 2D variables (horizontal cross sections)
	488
	489	REAL(wp) :: fill_value !< value for the _FillValue attribute
	490
	491	REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !<
	492	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which needs to be resorted for output
	493
	494	flag_nr = 0
	495	found = .TRUE.
	496	resorted = .FALSE.
	497	two_d = .FALSE.
	498
	499	SELECT CASE ( TRIM( variable ) )
	500
	501	CASE ( 'ti_xy', 'ti_xz', 'ti_yz' )
	502	IF ( av == 0 ) THEN
	503	to_be_resorted => ti
	504	ELSE
	505	IF ( .NOT. ALLOCATED( ti_av ) ) THEN
	506	ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
	507	ti_av = REAL( fill_value, KIND = wp )
	508	ENDIF
	509	to_be_resorted => ti_av
	510	ENDIF
	511	flag_nr = 0
	512
	513	IF ( mode == 'xy' ) grid = 'zu'
	514
	515	CASE ( 'uu_xy', 'uu_xz', 'uu_yz' )
	516	IF ( av == 0 ) THEN
	517	to_be_resorted => uu
	518	ELSE
	519	IF ( .NOT. ALLOCATED( uu_av ) ) THEN
	520	ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
	521	uu_av = REAL( fill_value, KIND = wp )
	522	ENDIF
	523	to_be_resorted => uu_av
	524	ENDIF
	525	flag_nr = 1
	526
	527	IF ( mode == 'xy' ) grid = 'zu'
	528
	529	CASE ( 'vv_xy', 'vv_xz', 'vv_yz' )
	530	IF ( av == 0 ) THEN
	531	to_be_resorted => vv
	532	ELSE
	533	IF ( .NOT. ALLOCATED( vv_av ) ) THEN
	534	ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
	535	vv_av = REAL( fill_value, KIND = wp )
	536	ENDIF
	537	to_be_resorted => vv_av
	538	ENDIF
	539	flag_nr = 2
	540
	541	IF ( mode == 'xy' ) grid = 'zu'
	542
	543	CASE ( 'ww_xy', 'ww_xz', 'ww_yz' )
	544	IF ( av == 0 ) THEN
	545	to_be_resorted => ww
	546	ELSE
	547	IF ( .NOT. ALLOCATED( ww_av ) ) THEN
	548	ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
	549	ww_av = REAL( fill_value, KIND = wp )
	550	ENDIF
	551	to_be_resorted => ww_av
	552	ENDIF
	553	flag_nr = 3
	554
	555	IF ( mode == 'xy' ) grid = 'zw'
	556
	557	CASE DEFAULT
	558	found = .FALSE.
	559	grid = 'none'
	560
	561	END SELECT
	562
	563	IF ( found .AND. .NOT. resorted ) THEN
	564	DO i = nxl, nxr
	565	DO j = nys, nyn
	566	DO k = nzb_do, nzt_do
	567	local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), &
	568	REAL( fill_value, KIND = wp ), &
	569	BTEST( wall_flags_0(k,j,i), flag_nr ) )
	570	ENDDO
	571	ENDDO
	572	ENDDO
[3994]	573	ENDIF
[4039]	574
	575	END SUBROUTINE doq_output_2d
	576
	577
	578	!------------------------------------------------------------------------------!
	579	!
	580	! Description:
	581	! ------------
	582	!> Subroutine defining 3D output variables
	583	!------------------------------------------------------------------------------!
	584	SUBROUTINE doq_output_3d( av, variable, found, local_pf, fill_value, nzb_do, &
	585	nzt_do )
	586
	587	IMPLICIT NONE
[3994]	588
[4039]	589	CHARACTER (LEN=*) :: variable !<
[3994]	590
[4039]	591	INTEGER(iwp) :: av !< index indicating averaged or instantaneous output
	592	INTEGER(iwp) :: flag_nr !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
	593	INTEGER(iwp) :: i !< index variable along x-direction
	594	INTEGER(iwp) :: j !< index variable along y-direction
	595	INTEGER(iwp) :: k !< index variable along z-direction
	596	INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0)
	597	INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d)
[3994]	598
[4039]	599	LOGICAL :: found !< true if variable is in list
	600	LOGICAL :: resorted !< true if array is resorted
[3994]	601
[4039]	602	REAL(wp) :: fill_value !< value for the _FillValue attribute
	603
	604	REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !<
	605	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which needs to be resorted for output
	606
	607	flag_nr = 0
	608	found = .TRUE.
	609	resorted = .FALSE.
	610
	611	SELECT CASE ( TRIM( variable ) )
	612
	613	CASE ( 'ti' )
	614	IF ( av == 0 ) THEN
	615	to_be_resorted => ti
	616	ELSE
	617	IF ( .NOT. ALLOCATED( ti_av ) ) THEN
	618	ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
	619	ti_av = REAL( fill_value, KIND = wp )
	620	ENDIF
	621	to_be_resorted => ti_av
	622	ENDIF
	623	flag_nr = 0
	624
	625	CASE ( 'uu' )
	626	IF ( av == 0 ) THEN
	627	to_be_resorted => uu
	628	ELSE
	629	IF ( .NOT. ALLOCATED( uu_av ) ) THEN
	630	ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
	631	uu_av = REAL( fill_value, KIND = wp )
	632	ENDIF
	633	to_be_resorted => uu_av
	634	ENDIF
	635	flag_nr = 1
	636
	637	CASE ( 'vv' )
	638	IF ( av == 0 ) THEN
	639	to_be_resorted => vv
	640	ELSE
	641	IF ( .NOT. ALLOCATED( vv_av ) ) THEN
	642	ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
	643	vv_av = REAL( fill_value, KIND = wp )
	644	ENDIF
	645	to_be_resorted => vv_av
	646	ENDIF
	647	flag_nr = 2
	648
	649	CASE ( 'ww' )
	650	IF ( av == 0 ) THEN
	651	to_be_resorted => ww
	652	ELSE
	653	IF ( .NOT. ALLOCATED( ww_av ) ) THEN
	654	ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
	655	ww_av = REAL( fill_value, KIND = wp )
	656	ENDIF
	657	to_be_resorted => ww_av
	658	ENDIF
	659	flag_nr = 3
	660
	661	CASE DEFAULT
	662	found = .FALSE.
	663
	664	END SELECT
	665
	666	IF ( found .AND. .NOT. resorted ) THEN
	667	DO i = nxl, nxr
	668	DO j = nys, nyn
	669	DO k = nzb_do, nzt_do
	670	local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), &
	671	REAL( fill_value, KIND = wp ), &
	672	BTEST( wall_flags_0(k,j,i), flag_nr ) )
	673	ENDDO
	674	ENDDO
	675	ENDDO
	676	ENDIF
	677
	678	END SUBROUTINE doq_output_3d
	679
[3994]	680	! Description:
	681	! ------------
[4039]	682	!> Resorts the user-defined output quantity with indices (k,j,i) to a
	683	!> temporary array with indices (i,j,k) for masked data output.
	684	!------------------------------------------------------------------------------!
[4069]	685	SUBROUTINE doq_output_mask( av, variable, found, local_pf, mid )
[4039]	686
	687	USE control_parameters
	688
	689	USE indices
	690
	691	USE surface_mod, &
	692	ONLY: get_topography_top_index_ji
	693
	694	IMPLICIT NONE
[3994]	695
[4039]	696	CHARACTER (LEN=*) :: variable !<
	697	CHARACTER (LEN=5) :: grid !< flag to distinquish between staggered grids
[3994]	698
[4039]	699	INTEGER(iwp) :: av !< index indicating averaged or instantaneous output
	700	INTEGER(iwp) :: flag_nr !< number of the topography flag (0: scalar, 1: u, 2: v, 3: w)
	701	INTEGER(iwp) :: i !< index variable along x-direction
	702	INTEGER(iwp) :: j !< index variable along y-direction
	703	INTEGER(iwp) :: k !< index variable along z-direction
[4069]	704	INTEGER(iwp) :: mid !< masked output running index
[4039]	705	INTEGER(iwp) :: topo_top_ind !< k index of highest horizontal surface
[3994]	706
[4039]	707	LOGICAL :: found !< true if variable is in list
	708	LOGICAL :: resorted !< true if array is resorted
[3994]	709
[4039]	710	REAL(wp), &
	711	DIMENSION(mask_size_l(mid,1),mask_size_l(mid,2),mask_size_l(mid,3)) :: &
	712	local_pf !<
	713	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to array which needs to be resorted for output
[3994]	714
[4039]	715	flag_nr = 0
	716	found = .TRUE.
	717	resorted = .FALSE.
	718	grid = 's'
	719
	720	SELECT CASE ( TRIM( variable ) )
	721
	722	CASE ( 'ti' )
	723	IF ( av == 0 ) THEN
	724	to_be_resorted => ti
	725	ELSE
	726	to_be_resorted => ti_av
	727	ENDIF
	728	grid = 's'
	729	flag_nr = 0
	730
	731	CASE ( 'uu' )
	732	IF ( av == 0 ) THEN
	733	to_be_resorted => uu
	734	ELSE
	735	to_be_resorted => uu_av
	736	ENDIF
	737	grid = 'u'
	738	flag_nr = 1
	739
	740	CASE ( 'vv' )
	741	IF ( av == 0 ) THEN
	742	to_be_resorted => vv
	743	ELSE
	744	to_be_resorted => vv_av
	745	ENDIF
	746	grid = 'v'
	747	flag_nr = 2
	748
	749	CASE ( 'ww' )
	750	IF ( av == 0 ) THEN
	751	to_be_resorted => ww
	752	ELSE
	753	to_be_resorted => ww_av
	754	ENDIF
	755	grid = 'w'
	756	flag_nr = 3
	757
	758
	759	CASE DEFAULT
	760	found = .FALSE.
	761
	762	END SELECT
	763
	764	IF ( found .AND. .NOT. resorted ) THEN
	765	IF ( .NOT. mask_surface(mid) ) THEN
	766	!
	767	!-- Default masked output
	768	DO i = 1, mask_size_l(mid,1)
	769	DO j = 1, mask_size_l(mid,2)
	770	DO k = 1, mask_size_l(mid,3)
	771	local_pf(i,j,k) = to_be_resorted(mask_k(mid,k), &
	772	mask_j(mid,j), &
	773	mask_i(mid,i))
	774	ENDDO
	775	ENDDO
	776	ENDDO
	777
	778	ELSE
	779	!
	780	!-- Terrain-following masked output
	781	DO i = 1, mask_size_l(mid,1)
	782	DO j = 1, mask_size_l(mid,2)
	783	!
	784	!-- Get k index of highest horizontal surface
	785	topo_top_ind = get_topography_top_index_ji( mask_j(mid,j), &
	786	mask_i(mid,i), &
	787	grid )
	788	!
	789	!-- Save output array
	790	DO k = 1, mask_size_l(mid,3)
	791	local_pf(i,j,k) = to_be_resorted( &
	792	MIN( topo_top_ind+mask_k(mid,k), &
	793	nzt+1 ), &
	794	mask_j(mid,j), &
	795	mask_i(mid,i) )
	796	ENDDO
	797	ENDDO
	798	ENDDO
	799
	800	ENDIF
	801	ENDIF
	802
	803	END SUBROUTINE doq_output_mask
	804
	805	!------------------------------------------------------------------------------!
	806	! Description:
	807	! ------------
	808	!> Allocate required arrays
	809	!------------------------------------------------------------------------------!
	810	SUBROUTINE doq_init
	811
[3994]	812	IMPLICIT NONE
[4039]	813
	814	INTEGER(iwp) :: ivar !< loop index over all 2d/3d/mask output quantities
	815	!
	816	!-- Next line is to avoid compiler warnings about unused variables
	817	IF ( timestep_number_at_prev_calc == 0 ) CONTINUE
[3994]	818
[4039]	819	initialized_diagnostic_output_quantities = .FALSE.
	820
	821	ivar = 1
	822
	823	DO WHILE ( ivar <= SIZE( do_all ) )
	824
	825	SELECT CASE ( TRIM( do_all(ivar) ) )
	826	!
	827	!-- Allocate array for 'turbulence intensity'
	828	CASE ( 'ti' )
	829	IF ( .NOT. ALLOCATED( ti ) ) THEN
	830	ALLOCATE( ti(nzb:nzt+1,nys:nyn,nxl:nxr) )
	831	ti = 0.0_wp
	832	ENDIF
	833	!
	834	!-- Allocate array for uu
	835	CASE ( 'uu' )
	836	IF ( .NOT. ALLOCATED( uu ) ) THEN
	837	ALLOCATE( uu(nzb:nzt+1,nys:nyn,nxl:nxr) )
	838	uu = 0.0_wp
	839	ENDIF
	840
	841	!
	842	!-- Allocate array for vv
	843	CASE ( 'vv' )
	844	IF ( .NOT. ALLOCATED( vv ) ) THEN
	845	ALLOCATE( vv(nzb:nzt+1,nys:nyn,nxl:nxr) )
	846	vv = 0.0_wp
	847	ENDIF
	848
	849	!
	850	!-- Allocate array for ww
	851	CASE ( 'ww' )
	852	IF ( .NOT. ALLOCATED( ww ) ) THEN
	853	ALLOCATE( ww(nzb:nzt+1,nys:nyn,nxl:nxr) )
	854	ww = 0.0_wp
	855	ENDIF
	856
	857	END SELECT
	858
	859	ivar = ivar + 1
	860	ENDDO
	861
	862	initialized_diagnostic_output_quantities = .TRUE.
	863
	864	END SUBROUTINE doq_init
	865
	866
	867	!--------------------------------------------------------------------------------------------------!
	868	! Description:
	869	! ------------
	870	!> Calculate diagnostic quantities
	871	!--------------------------------------------------------------------------------------------------!
	872	SUBROUTINE doq_calculate
	873
	874	IMPLICIT NONE
	875
[3994]	876	INTEGER(iwp) :: i, j, k !< grid loop index in x-, y-, z-direction
[4039]	877	INTEGER(iwp) :: ivar !< loop index over all 2d/3d/mask output quantities
[3994]	878
	879
	880	! CALL cpu_log( log_point(41), 'calculate_quantities', 'start' )
	881	!
	882	!-- Preparatory steps and initialization of output arrays
[4039]	883	IF ( .NOT. prepared_diagnostic_output_quantities ) CALL doq_prepare
	884	IF ( .NOT. initialized_diagnostic_output_quantities ) CALL doq_init
[3994]	885	!
[4039]	886	!-- Save timestep number to check in time_integration if doq_calculate
[3994]	887	!-- has been called already, since the CALL occurs at two locations, but the calculations need to be
	888	!-- done only once per timestep.
	889	timestep_number_at_prev_calc = current_timestep_number
	890
[4039]	891	ivar = 1
[3994]	892
[4039]	893	DO WHILE ( ivar <= SIZE( do_all ) )
[3994]	894
[4039]	895	SELECT CASE ( TRIM( do_all(ivar) ) )
[3994]	896	!
	897	!-- Calculate 'turbulence intensity' from rot[(u,v,w)] at scalar grid point
	898	CASE ( 'ti' )
	899	DO i = nxl, nxr
	900	DO j = nys, nyn
	901	DO k = nzb+1, nzt
[4039]	902	ti(k,j,i) = 0.25_wp * SQRT( &
	903	( ( w(k,j+1,i) + w(k-1,j+1,i) &
	904	- w(k,j-1,i) - w(k-1,j-1,i) ) * ddy &
	905	- ( v(k+1,j,i) + v(k+1,j+1,i) &
	906	- v(k-1,j,i) - v(k-1,j+1,i) ) * ddzu(k) )**2 &
	907	+ ( ( u(k+1,j,i) + u(k+1,j,i+1) &
	908	- u(k-1,j,i) - u(k-1,j,i+1) ) * ddzu(k) &
	909	- ( w(k,j,i+1) + w(k-1,j,i+1) &
	910	- w(k,j,i-1) - w(k-1,j,i-1) ) * ddx )**2 &
	911	+ ( ( v(k,j,i+1) + v(k,j+1,i+1) &
	912	- v(k,j,i-1) - v(k,j+1,i-1) ) * ddx &
	913	- ( u(k,j+1,i) + u(k,j+1,i+1) &
	914	- u(k,j-1,i) - u(k,j-1,i+1) ) * ddy )**2 ) &
	915	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0) )
[3994]	916	ENDDO
	917	ENDDO
[4039]	918	ENDDO
	919	!
	920	!-- uu
	921	CASE ( 'uu' )
	922	DO i = nxl, nxr
	923	DO j = nys, nyn
	924	DO k = nzb+1, nzt
	925	uu(k,j,i) = u(k,j,i) * u(k,j,i) &
	926	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 1) )
	927	ENDDO
	928	ENDDO
[3994]	929	ENDDO
[4039]	930	!
	931	!-- vv
	932	CASE ( 'vv' )
	933	DO i = nxl, nxr
	934	DO j = nys, nyn
	935	DO k = nzb+1, nzt
	936	vv(k,j,i) = v(k,j,i) * v(k,j,i) &
	937	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 2) )
	938	ENDDO
	939	ENDDO
	940	ENDDO
	941	!
	942	!-- ww
	943	CASE ( 'ww' )
	944	DO i = nxl, nxr
	945	DO j = nys, nyn
	946	DO k = nzb+1, nzt-1
	947	ww(k,j,i) = w(k,j,i) * w(k,j,i) &
	948	* MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 3) )
	949	ENDDO
	950	ENDDO
	951	ENDDO
[3994]	952
[4039]	953	END SELECT
[3994]	954
[4039]	955	ivar = ivar + 1
[3994]	956	ENDDO
	957
	958	! CALL cpu_log( log_point(41), 'calculate_quantities', 'stop' )
	959
[4039]	960	END SUBROUTINE doq_calculate
[3994]	961
	962
	963	!------------------------------------------------------------------------------!
	964	! Description:
	965	! ------------
	966	!> ...
	967	!------------------------------------------------------------------------------!
[4039]	968	SUBROUTINE doq_prepare
[3994]	969
	970
[4039]	971	USE control_parameters, &
[4069]	972	ONLY: do2d, do3d, domask, masks
[3994]	973
	974	IMPLICIT NONE
	975
	976	CHARACTER (LEN=varnamelength), DIMENSION(0:1,500) :: do2d_var = ' ' !<
	977	!< label array for 2d output quantities
	978
	979	INTEGER(iwp) :: av !< index defining type of output, av=0 instantaneous, av=1 averaged
	980	INTEGER(iwp) :: ivar !< loop index
	981	INTEGER(iwp) :: ivar_all !< loop index
	982	INTEGER(iwp) :: l !< index for cutting string
[4069]	983	INTEGER(iwp) :: mid !< masked output running index
[3994]	984
	985	prepared_diagnostic_output_quantities = .FALSE.
	986
	987	ivar = 1
	988	ivar_all = 1
	989
	990	DO av = 0, 1
	991	!
	992	!-- Remove _xy, _xz, or _yz from string
	993	l = MAX( 3, LEN_TRIM( do2d(av,ivar) ) )
[3998]	994	do2d_var(av,ivar)(1:l-3) = do2d(av,ivar)(1:l-3)
[3994]	995	!
[4039]	996	!-- Gather 2d output quantity names.
	997	!-- Check for double occurrence of output quantity, e.g. by _xy,
	998	!-- _yz, _xz.
[3994]	999	DO WHILE ( do2d_var(av,ivar)(1:1) /= ' ' )
	1000	IF ( .NOT. ANY( do_all == do2d_var(av,ivar) ) ) THEN
	1001	do_all(ivar_all) = do2d_var(av,ivar)
	1002	ENDIF
	1003	ivar = ivar + 1
	1004	ivar_all = ivar_all + 1
	1005	l = MAX( 3, LEN_TRIM( do2d(av,ivar) ) )
[3998]	1006	do2d_var(av,ivar)(1:l-3) = do2d(av,ivar)(1:l-3)
[3994]	1007	ENDDO
	1008
	1009	ivar = 1
	1010	!
[4039]	1011	!-- Gather 3d output quantity names
[3994]	1012	DO WHILE ( do3d(av,ivar)(1:1) /= ' ' )
	1013	do_all(ivar_all) = do3d(av,ivar)
	1014	ivar = ivar + 1
	1015	ivar_all = ivar_all + 1
	1016	ENDDO
	1017
	1018	ivar = 1
	1019	!
[4039]	1020	!-- Gather masked output quantity names. Also check for double output
	1021	!-- e.g. by different masks.
[3994]	1022	DO mid = 1, masks
	1023	DO WHILE ( domask(mid,av,ivar)(1:1) /= ' ' )
[4039]	1024	IF ( .NOT. ANY( do_all == domask(mid,av,ivar) ) ) THEN
	1025	do_all(ivar_all) = do2d_var(av,ivar)
	1026	ENDIF
[3994]	1027
	1028	ivar = ivar + 1
	1029	ivar_all = ivar_all + 1
	1030	ENDDO
	1031	ivar = 1
	1032	ENDDO
	1033
	1034	ENDDO
	1035
	1036	prepared_diagnostic_output_quantities = .TRUE.
	1037
[4039]	1038	END SUBROUTINE doq_prepare
	1039
	1040	!------------------------------------------------------------------------------!
	1041	! Description:
	1042	! ------------
	1043	!> Subroutine reads local (subdomain) restart data
	1044	!> Note: With the current structure reading of non-standard array is not
	1045	!> possible
	1046	!------------------------------------------------------------------------------!
	1047	! SUBROUTINE doq_rrd_local( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, &
	1048	! nxr_on_file, nynf, nync, nyn_on_file, nysf, &
	1049	! nysc, nys_on_file, tmp_3d_non_standard, found )
	1050	!
	1051	!
	1052	! USE control_parameters
	1053	!
	1054	! USE indices
	1055	!
	1056	! USE kinds
	1057	!
	1058	! USE pegrid
	1059	!
	1060	!
	1061	! IMPLICIT NONE
	1062	!
	1063	! INTEGER(iwp) :: k !<
	1064	! INTEGER(iwp) :: nxlc !<
	1065	! INTEGER(iwp) :: nxlf !<
	1066	! INTEGER(iwp) :: nxl_on_file !<
	1067	! INTEGER(iwp) :: nxrc !<
	1068	! INTEGER(iwp) :: nxrf !<
	1069	! INTEGER(iwp) :: nxr_on_file !<
	1070	! INTEGER(iwp) :: nync !<
	1071	! INTEGER(iwp) :: nynf !<
	1072	! INTEGER(iwp) :: nyn_on_file !<
	1073	! INTEGER(iwp) :: nysc !<
	1074	! INTEGER(iwp) :: nysf !<
	1075	! INTEGER(iwp) :: nys_on_file !<
	1076	!
	1077	! LOGICAL, INTENT(OUT) :: found
	1078	!
	1079	! REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: tmp_3d_non_standard !< temporary array for storing 3D data with non standard dimensions
	1080	! !
	1081	! !-- If temporary non-standard array for reading is already allocated,
	1082	! !-- deallocate it.
	1083	! IF ( ALLOCATED( tmp_3d_non_standard ) ) DEALLOCATE( tmp_3d_non_standard )
	1084	!
	1085	! found = .TRUE.
	1086	!
	1087	! SELECT CASE ( restart_string(1:length) )
	1088	!
	1089	! CASE ( 'ti_av' )
	1090	! IF ( .NOT. ALLOCATED( ti_av ) ) THEN
	1091	! ALLOCATE( ti_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
	1092	! ENDIF
	1093	! IF ( k == 1 ) THEN
	1094	! ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file, &
	1095	! nxl_on_file:nxr_on_file) )
	1096	! READ ( 13 ) tmp_3d_non_standard
	1097	! ENDIF
	1098	! ti_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
	1099	!
	1100	! CASE ( 'uu_av' )
	1101	! IF ( .NOT. ALLOCATED( uu_av ) ) THEN
	1102	! ALLOCATE( uu_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
	1103	! ENDIF
	1104	! IF ( k == 1 ) THEN
	1105	! ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file, &
	1106	! nxl_on_file:nxr_on_file) )
	1107	! READ ( 13 ) tmp_3d_non_standard
	1108	! ENDIF
	1109	! uu_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
	1110	!
	1111	! CASE ( 'vv_av' )
	1112	! IF ( .NOT. ALLOCATED( vv_av ) ) THEN
	1113	! ALLOCATE( vv_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
	1114	! ENDIF
	1115	! IF ( k == 1 ) THEN
	1116	! ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file, &
	1117	! nxl_on_file:nxr_on_file) )
	1118	! READ ( 13 ) tmp_3d_non_standard
	1119	! ENDIF
	1120	! vv_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
	1121	!
	1122	! CASE ( 'ww_av' )
	1123	! IF ( .NOT. ALLOCATED( ww_av ) ) THEN
	1124	! ALLOCATE( ww_av(nzb:nzt+1,nys:nyn,nxl:nxr) )
	1125	! ENDIF
	1126	! IF ( k == 1 ) THEN
	1127	! ALLOCATE( tmp_3d_non_standard(nzb:nzt+1,nys_on_file:nyn_on_file, &
	1128	! nxl_on_file:nxr_on_file) )
	1129	! READ ( 13 ) tmp_3d_non_standard
	1130	! ENDIF
	1131	! ww_av(:,nysc:nync,nxlc:nxrc) = tmp_3d_non_standard(:,nysf:nynf,nxlf:nxrf)
	1132	!
	1133	!
	1134	! CASE DEFAULT
	1135	!
	1136	! found = .FALSE.
	1137	!
	1138	! END SELECT
	1139	!
	1140	! END SUBROUTINE doq_rrd_local
	1141
	1142	!------------------------------------------------------------------------------!
	1143	! Description:
	1144	! ------------
	1145	!> Subroutine writes local (subdomain) restart data
	1146	!------------------------------------------------------------------------------!
	1147	SUBROUTINE doq_wrd_local
[3994]	1148
[4039]	1149
	1150	IMPLICIT NONE
	1151
	1152	IF ( ALLOCATED( ti_av ) ) THEN
	1153	CALL wrd_write_string( 'ti_av' )
	1154	WRITE ( 14 ) ti_av
	1155	ENDIF
	1156
	1157	IF ( ALLOCATED( uu_av ) ) THEN
	1158	CALL wrd_write_string( 'uu_av' )
	1159	WRITE ( 14 ) uu_av
	1160	ENDIF
	1161
	1162	IF ( ALLOCATED( vv_av ) ) THEN
	1163	CALL wrd_write_string( 'vv_av' )
	1164	WRITE ( 14 ) vv_av
	1165	ENDIF
	1166
	1167	IF ( ALLOCATED( ww_av ) ) THEN
	1168	CALL wrd_write_string( 'ww_av' )
	1169	WRITE ( 14 ) ww_av
	1170	ENDIF
	1171
	1172
	1173	END SUBROUTINE doq_wrd_local
	1174
	1175
	1176
[3994]	1177	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 |