Home

Context Navigation

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

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