Home

Context Navigation

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

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