Home

Context Navigation

source: palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90 @ 2742

Last change on this file since 2742 was 2716, checked in by kanani, 7 years ago
Correction of "Former revisions" section
Property svn:keywords set to `Id`
File size: 35.6 KB

Rev	Line
[2259]	1	!> @synthetic_turbulence_generator_mod.f90
	2	!------------------------------------------------------------------------------!
[2696]	3	! This file is part of the PALM model system.
[2259]	4	!
[2696]	5	! PALM is free software: you can redistribute it and/or modify it under the
[2259]	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	!
[2696]	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
[2259]	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 2017 Leibniz Universitaet Hannover
	18	!------------------------------------------------------------------------------!
	19	!
	20	! Current revisions:
	21	! -----------------
	22	!
	23	!
	24	! Former revisions:
	25	! -----------------
	26	! $Id: synthetic_turbulence_generator_mod.f90 2716 2017-12-29 16:35:59Z suehring $
[2716]	27	! Corrected "Former revisions" section
	28	!
	29	! 2696 2017-12-14 17:12:51Z kanani
	30	! Change in file header (GPL part)
	31	!
	32	! 2669 2017-12-06 16:03:27Z raasch
[2669]	33	! unit number for file containing turbulence generator data changed to 90
	34	! bugfix: preprocessor directives added for MPI specific code
	35	!
	36	! 2576 2017-10-24 13:49:46Z Giersch
[2576]	37	! Definition of a new function called stg_skip_var_list to skip module
	38	! parameters during reading restart data
	39	!
	40	! 2563 2017-10-19 15:36:10Z Giersch
[2563]	41	! stg_read_restart_data is called in stg_parin in the case of a restart run
	42	!
	43	! 2259 2017-06-08 09:09:11Z gronemeier
[2259]	44	! Initial revision
	45	!
	46	!
	47	!
	48	! Authors:
	49	! --------
	50	! @author Tobias Gronemeier, Atsushi Inagaki, Micha Gryschka, Christoph Knigge
	51	!
	52	!
	53	! Description:
	54	! ------------
	55	!> The module generates turbulence at the inflow boundary based on a method by
	56	!> Xie and Castro (2008) utilizing a Lund rotation (Lund, 1998) and a mass-flux
	57	!> correction by Kim et al. (2013).
	58	!> The turbulence is correlated based on length scales in y- and z-direction and
	59	!> a time scale for each velocity component. The profiles of length and time
	60	!> scales, mean u, v, w, e and pt, and all components of the Reynolds stress
	61	!> tensor are read from file STG_PROFILES.
	62	!>
	63	!> @todo test restart
	64	!> enable cyclic_fill
	65	!> implement turbulence generation for e and pt
	66	!> @note <Enter notes on the module>
	67	!> @bug Height information from input file is not used. Profiles from input
	68	!> must match with current PALM grid.
	69	!> Transformation of length scales to number of gridpoints does not
	70	!> consider grid stretching.
	71	!> In case of restart, velocity seeds differ from precursor run if a11,
	72	!> a22, or a33 are zero.
	73	!------------------------------------------------------------------------------!
	74	MODULE synthetic_turbulence_generator_mod
	75
	76
	77	USE arrays_3d, &
	78	ONLY: mean_inflow_profiles, u, v, w
	79
	80	USE constants, &
	81	ONLY: pi
	82
	83	USE control_parameters, &
	84	ONLY: initializing_actions, message_string, &
	85	synthetic_turbulence_generator
	86
	87	USE cpulog, &
	88	ONLY: cpu_log, log_point, log_point_s
	89
	90	USE indices, &
	91	ONLY: nbgp, nzb, nzt, nyng, nysg
	92
	93	USE kinds
	94
	95	USE MPI
	96
	97	USE pegrid, &
	98	ONLY: comm1dx, comm1dy, comm2d, id_inflow, ierr, myidx, pdims
	99
	100	USE transpose_indices, &
	101	ONLY: nzb_x, nzt_x
	102
	103
	104	IMPLICIT NONE
	105
	106	LOGICAL :: velocity_seed_initialized = .FALSE. !< true after first call of stg_main
	107	LOGICAL :: use_synthetic_turbulence_generator = .FALSE. !< switch to use synthetic turbulence generator
	108
	109	INTEGER(iwp) :: stg_type_yz !< MPI type for full z range
	110	INTEGER(iwp) :: stg_type_yz_small !< MPI type for small z range
	111	INTEGER(iwp) :: merg !< maximum length scale (in gp)
	112	INTEGER(iwp) :: mergp !< merg + nbgp
	113
	114	REAL(wp) :: mc_factor !< mass flux correction factor
	115
	116	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: displs !< displacement for MPI_GATHERV
	117	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: recv_count !< receive count for MPI_GATHERV
	118	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nuy !< length scale of u in y direction (in gp)
	119	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nuz !< length scale of u in z direction (in gp)
	120	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvy !< length scale of v in y direction (in gp)
	121	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvz !< length scale of v in z direction (in gp)
	122	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwy !< length scale of w in y direction (in gp)
	123	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwz !< length scale of w in z direction (in gp)
	124
	125	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: seed !< seed of random number for rn-generator
	126
	127	REAL(wp), DIMENSION(:), ALLOCATABLE :: a11 !< coefficient for Lund rotation
	128	REAL(wp), DIMENSION(:), ALLOCATABLE :: a21 !< coefficient for Lund rotation
	129	REAL(wp), DIMENSION(:), ALLOCATABLE :: a22 !< coefficient for Lund rotation
	130	REAL(wp), DIMENSION(:), ALLOCATABLE :: a31 !< coefficient for Lund rotation
	131	REAL(wp), DIMENSION(:), ALLOCATABLE :: a32 !< coefficient for Lund rotation
	132	REAL(wp), DIMENSION(:), ALLOCATABLE :: a33 !< coefficient for Lund rotation
	133	REAL(wp), DIMENSION(:), ALLOCATABLE :: tu !< Lagrangian time scale of u
	134	REAL(wp), DIMENSION(:), ALLOCATABLE :: tv !< Lagrangian time scale of v
	135	REAL(wp), DIMENSION(:), ALLOCATABLE :: tw !< Lagrangian time scale of w
	136
	137	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: buy !< filter function for u in y direction
	138	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: buz !< filter function for u in z direction
	139	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvy !< filter function for v in y direction
	140	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvz !< filter function for v in z direction
	141	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwy !< filter function for w in y direction
	142	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwz !< filter function for w in z direction
	143	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fu !< velocity seed for u
	144	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fuo !< velocity seed for u with new random number
	145	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fv !< velocity seed for v
	146	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fvo !< velocity seed for v with new random number
	147	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fw !< velocity seed for w
	148	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fwo !< velocity seed for w with new random number
	149
	150
	151	!
	152	!-- PALM interfaces:
	153	!-- Input parameter checks to be done in check_parameters
	154	INTERFACE stg_check_parameters
	155	MODULE PROCEDURE stg_check_parameters
	156	END INTERFACE stg_check_parameters
	157
	158	!
	159	!-- Calculate filter functions
	160	INTERFACE stg_filter_func
	161	MODULE PROCEDURE stg_filter_func
	162	END INTERFACE stg_filter_func
	163
	164	!
	165	!-- Generate velocity seeds
	166	INTERFACE stg_generate_seed_yz
	167	MODULE PROCEDURE stg_generate_seed_yz
	168	END INTERFACE stg_generate_seed_yz
	169
	170	!
	171	!-- Output of information to the header file
	172	INTERFACE stg_header
	173	MODULE PROCEDURE stg_header
	174	END INTERFACE stg_header
	175
	176	!
	177	!-- Initialization actions
	178	INTERFACE stg_init
	179	MODULE PROCEDURE stg_init
	180	END INTERFACE stg_init
	181
	182	!
	183	!-- Main procedure of synth. turb. gen.
	184	INTERFACE stg_main
	185	MODULE PROCEDURE stg_main
	186	END INTERFACE stg_main
	187
	188	!
	189	!-- Reading of NAMELIST parameters
	190	INTERFACE stg_parin
	191	MODULE PROCEDURE stg_parin
	192	END INTERFACE stg_parin
	193
	194	!
[2576]	195	!-- Skipping of parameters for restart runs
	196	INTERFACE stg_skip_var_list
	197	MODULE PROCEDURE stg_skip_var_list
	198	END INTERFACE stg_skip_var_list
	199
	200	!
[2259]	201	!-- Reading of parameters for restart runs
	202	INTERFACE stg_read_restart_data
	203	MODULE PROCEDURE stg_read_restart_data
	204	END INTERFACE stg_read_restart_data
	205
	206	!
	207	!-- Writing of binary output for restart runs
	208	INTERFACE stg_write_restart_data
	209	MODULE PROCEDURE stg_write_restart_data
	210	END INTERFACE stg_write_restart_data
	211
	212	SAVE
	213
	214	PRIVATE
	215
	216	!
	217	!-- Public interfaces
	218	PUBLIC stg_check_parameters, stg_header, stg_init, stg_main, stg_parin, &
[2576]	219	stg_write_restart_data, stg_skip_var_list
[2259]	220
	221	!
	222	!-- Public variables
	223	PUBLIC use_synthetic_turbulence_generator
	224
	225
	226	CONTAINS
	227
	228
	229	!------------------------------------------------------------------------------!
	230	! Description:
	231	! ------------
	232	!> Check parameters routine for synthetic turbulence generator
	233	!------------------------------------------------------------------------------!
	234	SUBROUTINE stg_check_parameters
	235
	236
	237	USE control_parameters, &
	238	ONLY: bc_lr, bc_ns, turbulent_inflow
	239
	240
	241	IMPLICIT NONE
	242
	243	IF ( use_synthetic_turbulence_generator ) THEN
	244
	245	IF ( INDEX( initializing_actions, 'set_constant_profiles' ) == 0 .AND. &
	246	INDEX( initializing_actions, 'read_restart_data' ) == 0 ) THEN
	247	message_string = 'Using synthetic turbulence generator ' // &
	248	'requires &initializing_actions = ' // &
	249	'"set_constant_profiles" or "read_restart_data"'
	250	CALL message( 'stg_check_parameters', 'PA0015', 1, 2, 0, 6, 0 )
	251	ENDIF
	252
	253	IF ( bc_lr /= 'dirichlet/radiation' ) THEN
	254	message_string = 'Using synthetic turbulence generator ' // &
	255	'requires &bc_lr = "dirichlet/radiation"'
	256	CALL message( 'stg_check_parameters', 'PA0035', 1, 2, 0, 6, 0 )
	257	ENDIF
	258	IF ( bc_ns /= 'cyclic' ) THEN
	259	message_string = 'Using synthetic turbulence generator ' // &
	260	'requires &bc_ns = "cyclic"'
	261	CALL message( 'stg_check_parameters', 'PA0037', 1, 2, 0, 6, 0 )
	262	ENDIF
	263	IF ( turbulent_inflow ) THEN
	264	message_string = 'Using synthetic turbulence generator ' // &
	265	'in combination &with turbulent_inflow = .T. ' // &
	266	'is not allowed'
	267	CALL message( 'stg_check_parameters', 'PA0039', 1, 2, 0, 6, 0 )
	268	ENDIF
	269
	270	ENDIF
	271
	272	END SUBROUTINE stg_check_parameters
	273
	274
	275	!------------------------------------------------------------------------------!
	276	! Description:
	277	! ------------
	278	!> Header output for synthetic turbulence generator
	279	!------------------------------------------------------------------------------!
	280	SUBROUTINE stg_header ( io )
	281
	282
	283	IMPLICIT NONE
	284
	285	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
	286
	287	!
	288	!-- Write synthetic turbulence generator Header
	289	WRITE( io, 1 )
	290	IF ( use_synthetic_turbulence_generator ) THEN
	291	WRITE( io, 2 )
	292	ELSE
	293	WRITE( io, 3 )
	294	ENDIF
	295
	296	1 FORMAT (//' Synthetic turbulence generator information:'/ &
	297	' ------------------------------------------'/)
	298	2 FORMAT (' synthetic turbulence generator switched on')
	299	3 FORMAT (' synthetic turbulence generator switched off')
	300
	301	END SUBROUTINE stg_header
	302
	303
	304	!------------------------------------------------------------------------------!
	305	! Description:
	306	! ------------
	307	!> Initialization of the synthetic turbulence generator
	308	!------------------------------------------------------------------------------!
	309	SUBROUTINE stg_init
	310
	311
	312	USE arrays_3d, &
	313	ONLY: u_init, v_init
	314
	315	USE control_parameters, &
	316	ONLY: coupling_char, dz, e_init
	317
	318	USE grid_variables, &
	319	ONLY: ddy
	320
	321
	322	IMPLICIT NONE
	323
[2669]	324	#if defined( __parallel )
[2259]	325	INTEGER(KIND=MPI_ADDRESS_KIND) :: extent !< extent of new MPI type
	326	INTEGER(KIND=MPI_ADDRESS_KIND) :: tob=0 !< dummy variable
[2669]	327	#endif
[2259]	328
	329	INTEGER(iwp) :: j !> loop index
	330	INTEGER(iwp) :: k !< index
	331	INTEGER(iwp) :: newtype !< dummy MPI type
	332	INTEGER(iwp) :: realsize !< size of REAL variables
	333	INTEGER(iwp) :: nseed !< dimension of random number seed
	334	INTEGER(iwp) :: startseed = 1234567890 !< start seed for random number generator
	335
	336	!
	337	!-- Dummy variables used for reading profiles
	338	REAL(wp) :: d1 !< u profile
	339	REAL(wp) :: d2 !< v profile
	340	REAL(wp) :: d3 !< w profile
	341	REAL(wp) :: d5 !< e profile
	342	REAL(wp) :: d11 !< vertical interpolation for a11
	343	REAL(wp) :: d21 !< vertical interpolation for a21
	344	REAL(wp) :: d22 !< vertical interpolation for a22
	345	REAL(wp) :: luy !< length scale for u in y direction
	346	REAL(wp) :: luz !< length scale for u in z direction
	347	REAL(wp) :: lvy !< length scale for v in y direction
	348	REAL(wp) :: lvz !< length scale for v in z direction
	349	REAL(wp) :: lwy !< length scale for w in y direction
	350	REAL(wp) :: lwz !< length scale for w in z direction
	351	REAL(wp) :: zz !< height
	352
	353	REAL(wp),DIMENSION(nzb:nzt+1) :: r11 !< Reynolds parameter
	354	REAL(wp),DIMENSION(nzb:nzt+1) :: r21 !< Reynolds parameter
	355	REAL(wp),DIMENSION(nzb:nzt+1) :: r22 !< Reynolds parameter
	356	REAL(wp),DIMENSION(nzb:nzt+1) :: r31 !< Reynolds parameter
	357	REAL(wp),DIMENSION(nzb:nzt+1) :: r32 !< Reynolds parameter
	358	REAL(wp),DIMENSION(nzb:nzt+1) :: r33 !< Reynolds parameter
	359
	360
	361	#if defined( __parallel )
	362	CALL MPI_BARRIER( comm2d, ierr )
	363	#endif
	364
	365	CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
	366
	367	IF ( .NOT. ALLOCATED( mean_inflow_profiles ) ) &
	368	ALLOCATE( mean_inflow_profiles(nzb:nzt+1,5) )
	369
	370	ALLOCATE ( a11(nzb:nzt+1), a21(nzb:nzt+1), a22(nzb:nzt+1), &
	371	a31(nzb:nzt+1), a32(nzb:nzt+1), a33(nzb:nzt+1), &
	372	nuy(nzb:nzt+1), nuz(nzb:nzt+1), nvy(nzb:nzt+1), &
	373	nvz(nzb:nzt+1), nwy(nzb:nzt+1), nwz(nzb:nzt+1), &
	374	tu(nzb:nzt+1), tv(nzb:nzt+1), tw(nzb:nzt+1) )
	375
[2669]	376	#if defined( __parallel )
[2259]	377	!
	378	!-- Define MPI type used in stg_generate_seed_yz to gather vertical splitted
	379	!-- velocity seeds
	380	CALL MPI_TYPE_SIZE( MPI_REAL, realsize, ierr )
	381	extent = 1 * realsize
	382
	383	! stg_type_yz: yz-slice with vertical bounds nzb:nzt+1
	384	CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt-nzb+2,nyng-nysg+1], &
	385	[1,nyng-nysg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
	386	CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_yz, ierr )
	387	CALL MPI_TYPE_COMMIT( stg_type_yz, ierr )
	388	CALL MPI_TYPE_FREE( newtype, ierr )
	389
	390	! stg_type_yz_small: yz-slice with vertical bounds nzb_x:nzt_x+1
	391	CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt_x-nzb_x+2,nyng-nysg+1], &
	392	[1,nyng-nysg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
	393	CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_yz_small, ierr )
	394	CALL MPI_TYPE_COMMIT( stg_type_yz_small, ierr )
	395	CALL MPI_TYPE_FREE( newtype, ierr )
	396
	397	! receive count and displacement for MPI_GATHERV in stg_generate_seed_yz
	398	ALLOCATE( recv_count(pdims(1)), displs(pdims(1)) )
	399
	400	recv_count = nzt_x-nzb_x + 1
	401	recv_count(pdims(1)) = recv_count(pdims(1)) + 1
	402
	403	DO j = 1, pdims(1)
	404	displs(j) = 0 + (nzt_x-nzb_x+1) * (j-1)
	405	ENDDO
[2669]	406	#endif
[2259]	407
	408	!
	409	!-- Define seed of random number
	410	CALL RANDOM_SEED()
	411	CALL RANDOM_SEED( size=nseed )
	412	ALLOCATE( seed(1:nseed) )
	413	DO j = 1, nseed
	414	seed(j) = startseed + j
	415	ENDDO
	416	CALL RANDOM_SEED(put=seed)
	417
	418	!
	419	!-- Read inflow profile
	420	!-- Height levels of profiles in input profile is as follows:
	421	!-- zu: luy, luz, tu, lvy, lvz, tv, r11, r21, r22, d1, d2, d5
	422	!-- zw: lwy, lwz, tw, r31, r32, r33, d3
	423	!-- WARNING: zz is not used at the moment
[2669]	424	OPEN( 90, FILE='STG_PROFILES'//TRIM( coupling_char ), STATUS='OLD', &
[2259]	425	FORM='FORMATTED')
	426
	427	! Skip header
[2669]	428	READ( 90, * )
[2259]	429
	430	DO k = nzb, nzt+1
[2669]	431	READ( 90, * ) zz, luy, luz, tu(k), lvy, lvz, tv(k), lwy, lwz, tw(k), &
[2259]	432	r11(k), r21(k), r22(k), r31(k), r32(k), r33(k), &
	433	d1, d2, d3, d5
	434
	435	!
	436	!-- Convert length scales from meter to number of grid points
	437	nuy(k) = INT( luy * ddy )
	438	nuz(k) = INT( luz / dz )
	439	nvy(k) = INT( lvy * ddy )
	440	nvz(k) = INT( lvz / dz )
	441	nwy(k) = INT( lwy * ddy )
	442	nwz(k) = INT( lwz / dz )
	443
	444	!
	445	!-- Save Mean inflow profiles
	446	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
	447	mean_inflow_profiles(k,1) = d1
	448	mean_inflow_profiles(k,2) = d2
	449	! mean_inflow_profiles(k,4) = d4
	450	mean_inflow_profiles(k,5) = d5
	451	ENDIF
	452	ENDDO
	453
[2669]	454	CLOSE( 90 )
[2259]	455
	456	!
	457	!-- Assign initial profiles
	458	u_init = mean_inflow_profiles(:,1)
	459	v_init = mean_inflow_profiles(:,2)
	460	! pt_init = mean_inflow_profiles(:,4)
	461	e_init = MAXVAL( mean_inflow_profiles(:,5) )
	462
	463	!
	464	!-- Calculate coefficient matrix from Reynolds stress tensor (Lund rotation)
	465	DO k = nzb, nzt+1
	466	IF ( r11(k) > 0.0_wp ) THEN
	467	a11(k) = SQRT( r11(k) )
	468	a21(k) = r21(k) / a11(k)
	469	ELSE
	470	a11(k) = 0.0_wp
	471	a21(k) = 0.0_wp
	472	ENDIF
	473
	474	a22(k) = r22(k) - a21(k)**2
	475	IF ( a22(k) > 0.0_wp ) THEN
	476	a22(k) = SQRT( a22(k) )
	477	a32(k) = ( r32(k) - a21(k) * a31(k) ) / a22(k)
	478	ELSE
	479	a22(k) = 0.0_wp
	480	a32(k) = 0.0_wp
	481	ENDIF
	482
	483	!
	484	!-- a31, a32, a33 must be calculated with interpolated a11, a21, a22 (d11,
	485	!-- d21, d22) because of different vertical grid
	486	IF ( k .le. nzt ) THEN
	487	d11 = 0.5_wp * ( r11(k) + r11(k+1) )
	488	IF ( d11 > 0.0_wp ) THEN
	489	d11 = SQRT( d11 )
	490	d21 = ( 0.5_wp * ( r21(k) + r21(k+1) ) ) / d11
	491	a31(k) = r31(k) / d11
	492	ELSE
	493	d21 = 0.0_wp
	494	a31(k) = 0.0_wp
	495	ENDIF
	496
	497	d22 = 0.5_wp * ( r22(k) + r22(k+1) ) - d21 ** 2
	498	IF ( d22 > 0.0_wp ) THEN
	499	a32(k) = ( r32(k) - d21 * a31(k) ) / SQRT( d22 )
	500	ELSE
	501	a32(k) = 0.0_wp
	502	ENDIF
	503
	504	a33(k) = r33(k) - a31(k) 2 - a32(k) 2
	505	IF ( a33(k) > 0.0_wp ) THEN
	506	a33(k) = SQRT( a33(k) )
	507	ELSE
	508	a33(k) = 0.0_wp
	509	ENDIF
	510	ELSE
	511	a31(k) = a31(k-1)
	512	a32(k) = a32(k-1)
	513	a33(k) = a33(k-1)
	514	ENDIF
	515
	516	ENDDO
	517
	518	!
	519	!-- Define the size of the filter functions and allocate them.
	520	merg = 0
	521
	522	! arrays must be large enough to cover the largest length scale
	523	DO k = nzb, nzt+1
	524	j = MAX( ABS(nuy(k)), ABS(nuz(k)), &
	525	ABS(nvy(k)), ABS(nvz(k)), &
	526	ABS(nwy(k)), ABS(nwz(k)) )
	527	IF ( j > merg ) merg = j
	528	ENDDO
	529
	530	merg = 2 * merg
	531	mergp = merg + nbgp
	532
	533	ALLOCATE ( buy(-merg:merg,nzb:nzt+1), buz(-merg:merg,nzb:nzt+1), &
	534	bvy(-merg:merg,nzb:nzt+1), bvz(-merg:merg,nzb:nzt+1), &
	535	bwy(-merg:merg,nzb:nzt+1), bwz(-merg:merg,nzb:nzt+1) )
	536
	537	!
	538	!-- Allocate velocity seeds
	539	ALLOCATE ( fu( nzb:nzt+1,nysg:nyng), fuo(nzb:nzt+1,nysg:nyng), &
	540	fv( nzb:nzt+1,nysg:nyng), fvo(nzb:nzt+1,nysg:nyng), &
	541	fw( nzb:nzt+1,nysg:nyng), fwo(nzb:nzt+1,nysg:nyng) )
	542
	543	fu = 0._wp
	544	fuo = 0._wp
	545	fv = 0._wp
	546	fvo = 0._wp
	547	fw = 0._wp
	548	fwo = 0._wp
	549
	550	!
	551	!-- Create filter functions
	552	CALL stg_filter_func( nuy, buy ) !filter uy
	553	CALL stg_filter_func( nuz, buz ) !filter uz
	554	CALL stg_filter_func( nvy, bvy ) !filter vy
	555	CALL stg_filter_func( nvz, bvz ) !filter vz
	556	CALL stg_filter_func( nwy, bwy ) !filter wy
	557	CALL stg_filter_func( nwz, bwz ) !filter wz
	558
	559	#if defined( __parallel )
	560	CALL MPI_BARRIER( comm2d, ierr )
	561	#endif
	562
	563	!
	564	!-- In case of restart, calculate velocity seeds fu, fv, fw from former
	565	! time step.
	566	! Bug: fu, fv, fw are different in those heights where a11, a22, a33
	567	! are 0 compared to the prerun. This is mostly for k=nzt+1.
	568	IF ( TRIM( initializing_actions ) == 'read_restart_data' ) THEN
	569	IF ( myidx == id_inflow ) THEN
	570	DO j = nysg, nyng
	571	DO k = nzb, nzt+1
	572
	573	IF ( a11(k) .NE. 0._wp ) THEN
	574	fu(k,j) = ( u(k,j,-1) / mc_factor - &
	575	mean_inflow_profiles(k,1) ) / a11(k)
	576	ELSE
	577	fu(k,j) = 0._wp
	578	ENDIF
	579
	580	IF ( a22(k) .NE. 0._wp ) THEN
	581	fv(k,j) = ( v(k,j,-1) / mc_factor - a21(k) * fu(k,j) - &
	582	mean_inflow_profiles(k,2) ) / a22(k)
	583	ELSE
	584	fv(k,j) = 0._wp
	585	ENDIF
	586
	587	IF ( a33(k) .NE. 0._wp ) THEN
	588	fw(k,j) = ( w(k,j,-1) / mc_factor - a31(k) * fu(k,j) - &
	589	a32(k) * fv(k,j) ) / a33(k)
	590	ELSE
	591	fw = 0._wp
	592	ENDIF
	593
	594	ENDDO
	595	ENDDO
	596	ENDIF
	597	ENDIF
	598
	599	CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
	600
	601	END SUBROUTINE stg_init
	602
	603
	604	!------------------------------------------------------------------------------!
	605	! Description:
	606	! ------------
	607	!> Calculate filter function bxx from length scale nxx following Eg.9 and 10
	608	!> (Xie and Castro, 2008)
	609	!------------------------------------------------------------------------------!
	610	SUBROUTINE stg_filter_func( nxx, bxx )
	611
	612
	613	IMPLICIT NONE
	614
	615	INTEGER(iwp) :: k !< loop index
	616	INTEGER(iwp) :: n_k !< length scale nXX in height k
	617	INTEGER(iwp) :: n_k2 !< n_k * 2
	618	INTEGER(iwp) :: nf !< index for length scales
	619
	620	REAL(wp) :: bdenom !< denominator for filter functions bXX
	621	REAL(wp) :: qsi = 1.0_wp !< minimization factor
	622
	623	INTEGER(iwp), DIMENSION(:) :: nxx(nzb:nzt+1) !< length scale (in gp)
	624
	625	REAL(wp), DIMENSION(:,:) :: bxx(-merg:merg,nzb:nzt+1) !< filter function
	626
	627
	628	bxx = 0.0_wp
	629
	630	DO k = nzb, nzt+1
	631	bdenom = 0.0_wp
	632	n_k = nxx(k)
	633	IF ( n_k /= 0 ) THEN
	634	n_k2 = n_k * 2
	635
	636	!
	637	!-- ( Eq.10 )^2
	638	DO nf = -n_k2, n_k2
	639	bdenom = bdenom + EXP( -qsi * pi * ABS(nf) / n_k )**2
	640	ENDDO
	641
	642	!
	643	!-- ( Eq.9 )
	644	bdenom = SQRT( bdenom )
	645	DO nf = -n_k2, n_k2
	646	bxx(nf,k) = EXP( -qsi * pi * ABS(nf) / n_k ) / bdenom
	647	ENDDO
	648	ENDIF
	649	ENDDO
	650
	651	END SUBROUTINE stg_filter_func
	652
	653
	654	!------------------------------------------------------------------------------!
	655	! Description:
	656	! ------------
	657	!> Parin for &stg_par for synthetic turbulence generator
	658	!------------------------------------------------------------------------------!
	659	SUBROUTINE stg_parin
	660
	661
	662	IMPLICIT NONE
	663
	664	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
	665
	666
	667	NAMELIST /stg_par/ use_synthetic_turbulence_generator
	668
	669	line = ' '
	670
	671	!
	672	!-- Try to find stg package
	673	REWIND ( 11 )
	674	line = ' '
	675	DO WHILE ( INDEX( line, '&stg_par' ) == 0 )
	676	READ ( 11, '(A)', END=10 ) line
	677	ENDDO
	678	BACKSPACE ( 11 )
	679
	680	!
	681	!-- Read namelist
	682	READ ( 11, stg_par )
	683
	684	!
	685	!-- Set flag that indicates that the synthetic turbulence generator is switched
	686	!-- on
	687	synthetic_turbulence_generator = .TRUE.
	688
[2563]	689	IF ( TRIM( initializing_actions ) == 'read_restart_data' ) THEN
	690	CALL stg_read_restart_data
	691	ENDIF
	692
[2259]	693	10 CONTINUE
	694
	695	END SUBROUTINE stg_parin
	696
	697
	698	!------------------------------------------------------------------------------!
	699	! Description:
	700	! ------------
[2576]	701	!> Skipping the stg variables from restart-file (binary format).
	702	!------------------------------------------------------------------------------!
	703	SUBROUTINE stg_skip_var_list
	704
	705	IMPLICIT NONE
	706
	707	CHARACTER (LEN=1) :: cdum
	708	CHARACTER (LEN=30) :: variable_chr
	709
	710	READ ( 13 ) variable_chr
	711
	712	DO WHILE ( TRIM( variable_chr ) /= '* end stg module *' )
	713
	714	READ ( 13 ) cdum
	715	READ ( 13 ) variable_chr
	716
	717	ENDDO
	718
	719	END SUBROUTINE stg_skip_var_list
	720
	721
	722	!------------------------------------------------------------------------------!
	723	! Description:
	724	! ------------
[2259]	725	!> This routine reads the respective restart data.
	726	!------------------------------------------------------------------------------!
	727	SUBROUTINE stg_read_restart_data
	728
	729
	730	IMPLICIT NONE
	731
	732	CHARACTER (LEN=30) :: variable_chr !< dummy variable to read string
	733
	734
	735	READ ( 13 ) variable_chr
	736	DO WHILE ( TRIM( variable_chr ) /= '* end stg module *' )
	737
	738	SELECT CASE ( TRIM( variable_chr ) )
	739
	740	CASE ( 'use_synthetic_turbulence_generator ' )
	741	READ ( 13 ) use_synthetic_turbulence_generator
	742	CASE ( 'mc_factor' )
	743	READ ( 13 ) mc_factor
	744
	745	END SELECT
	746
	747	READ ( 13 ) variable_chr
	748
	749	ENDDO
	750
	751	END SUBROUTINE stg_read_restart_data
	752
	753
	754	!------------------------------------------------------------------------------!
	755	! Description:
	756	! ------------
	757	!> This routine writes the respective restart data.
	758	!------------------------------------------------------------------------------!
	759	SUBROUTINE stg_write_restart_data
	760
	761
	762	IMPLICIT NONE
	763
	764	WRITE ( 14 ) 'use_synthetic_turbulence_generator '
	765	WRITE ( 14 ) use_synthetic_turbulence_generator
	766	WRITE ( 14 ) 'mc_factor '
	767	WRITE ( 14 ) mc_factor
	768
	769	WRITE ( 14 ) '* end stg module * '
	770
	771	END SUBROUTINE stg_write_restart_data
	772
	773
	774	!------------------------------------------------------------------------------!
	775	! Description:
	776	! ------------
	777	!> Create turbulent inflow fields for u, v, w with prescribed length scales and
	778	!> Reynolds stress tensor after a method of Xie and Castro (2008), modified
	779	!> following suggestions of Kim et al. (2013), and using a Lund rotation
	780	!> (Lund, 1998).
	781	!------------------------------------------------------------------------------!
	782	SUBROUTINE stg_main
	783
	784
	785	USE arrays_3d, &
	786	ONLY: dzw
	787
	788	USE control_parameters, &
	789	ONLY: dt_3d, intermediate_timestep_count, simulated_time, &
	790	volume_flow_initial
	791
	792	USE indices, &
	793	ONLY: nyn, nys
	794
	795	USE statistics, &
	796	ONLY: weight_substep
	797
	798
	799	IMPLICIT NONE
	800
	801	INTEGER(iwp) :: j !< loop index in y-direction
	802	INTEGER(iwp) :: k !< loop index in z-direction
	803
	804	REAL(wp) :: dt_stg !< wheighted subtimestep
	805	REAL(wp) :: mc_factor_l !< local mass flux correction factor
	806
	807	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,5,nbgp) :: inflow_dist !< disturbance for inflow profiles
	808
	809
	810	CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
	811
	812	!
	813	!-- Calculate time step which is needed for filter functions
	814	dt_stg = dt_3d * weight_substep(intermediate_timestep_count)
	815
	816	!
	817	!-- Initial value of fu, fv, fw
	818	IF ( simulated_time == 0.0_wp .AND. .NOT. velocity_seed_initialized ) THEN
	819	CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu )
	820	CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv )
	821	CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw )
	822	velocity_seed_initialized = .TRUE.
	823	ENDIF
	824	!
	825	!-- New set of fu, fv, fw
	826	CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo )
	827	CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo )
	828	CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo )
	829
	830	IF ( myidx == id_inflow ) THEN
	831
	832	DO j = nysg, nyng
	833	DO k = nzb, nzt + 1
	834	!
	835	!-- Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
	836	IF ( tu(k) == 0.0_wp ) THEN
	837	fu(k,j) = fuo(k,j)
	838	ELSE
	839	fu(k,j) = fu(k,j) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) + &
	840	fuo(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )
	841	ENDIF
	842
	843	IF ( tv(k) == 0.0_wp ) THEN
	844	fv(k,j) = fvo(k,j)
	845	ELSE
	846	fv(k,j) = fv(k,j) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) + &
	847	fvo(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
	848	ENDIF
	849
	850	IF ( tw(k) == 0.0_wp ) THEN
	851	fw(k,j) = fwo(k,j)
	852	ELSE
	853	fw(k,j) = fw(k,j) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) + &
	854	fwo(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
	855	ENDIF
	856	!
	857	!-- Lund rotation following Eq. 17 in Xie and Castro (2008).
	858	!-- Additional factors are added to improve the variance of v and w
	859	IF( k == 0 ) THEN
	860	inflow_dist(k,j,1,1) = 0.0_wp
	861	inflow_dist(k,j,2,1) = 0.0_wp
	862	inflow_dist(k,j,3,1) = 0.0_wp
	863	inflow_dist(k,j,4,1) = 0.0_wp
	864	inflow_dist(k,j,5,1) = 0.0_wp
	865	ELSE
	866	inflow_dist(k,j,1,1) = a11(k) * fu(k,j)
	867	!experimental test of 1.2
	868	inflow_dist(k,j,2,1) = ( SQRT( a22(k) / MAXVAL(a22) ) &
	869	* 1.2_wp ) &
	870	* ( a21(k) * fu(k,j) &
	871	+ a22(k) * fv(k,j) )
	872	inflow_dist(k,j,3,1) = ( SQRT(a33(k) / MAXVAL(a33) ) &
	873	* 1.3_wp ) &
	874	* ( a31(k) * fu(k,j) &
	875	+ a32(k) * fv(k,j) &
	876	+ a33(k) * fw(k,j) )
	877	! Calculation for pt and e not yet implemented
	878	inflow_dist(k,j,4,1) = 0.0_wp
	879	inflow_dist(k,j,5,1) = 0.0_wp
	880	ENDIF
	881
	882	ENDDO
	883	ENDDO
	884
	885	!
	886	!-- Mass flux correction following Kim et al. (2013)
	887	!-- This correction factor insures that the mass flux is preserved at the
	888	!-- inflow boundary
	889	mc_factor_l = 0.0_wp
	890	mc_factor = 0.0_wp
	891	DO j = nys, nyn
	892	DO k = nzb+1, nzt
	893	mc_factor_l = mc_factor_l + dzw(k) * &
	894	( mean_inflow_profiles(k,1) + inflow_dist(k,j,1,1) )
	895	ENDDO
	896	ENDDO
	897
	898	#if defined( __parallel )
	899	CALL MPI_ALLREDUCE( mc_factor_l, mc_factor, &
	900	1, MPI_REAL, MPI_SUM, comm1dy, ierr )
	901	#else
	902	mc_factor = mc_factor_l
	903	#endif
	904
	905	mc_factor = volume_flow_initial(1) / mc_factor
	906
	907	!
	908	!-- Add disturbance at the inflow
	909	DO j = nysg, nyng
	910	DO k = nzb, nzt+1
	911	u(k,j,-nbgp+1:0) = ( mean_inflow_profiles(k,1) + &
	912	inflow_dist(k,j,1,1) ) * mc_factor
	913	v(k,j,-nbgp:-1) = ( mean_inflow_profiles(k,2) + &
	914	inflow_dist(k,j,2,1) ) * mc_factor
	915	w(k,j,-nbgp:-1) = inflow_dist(k,j,3,1) * mc_factor
	916	ENDDO
	917	ENDDO
	918
	919	ENDIF
	920
	921	CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
	922
	923	END SUBROUTINE stg_main
	924
	925
	926	!------------------------------------------------------------------------------!
	927	! Description:
	928	! ------------
	929	!> Generate a set of random number rand_it wich is equal on each PE
	930	!> and calculate the velocity seed f_n.
	931	!> f_n is splitted in vertical direction by the number of PEs in x-direction and
	932	!> and each PE calculates a vertical subsection of f_n. At the the end, all
	933	!> parts are collected to form the full array.
	934	!------------------------------------------------------------------------------!
	935	SUBROUTINE stg_generate_seed_yz( n_y, n_z, b_y, b_z, f_n )
	936
	937
	938	USE indices, &
	939	ONLY: ny
	940
	941
	942	IMPLICIT NONE
	943
	944	INTEGER(iwp) :: j !< loop index in y-direction
	945	INTEGER(iwp) :: jj !< loop index in y-direction
	946	INTEGER(iwp) :: k !< loop index in z-direction
	947	INTEGER(iwp) :: kk !< loop index in z-direction
	948	INTEGER(iwp) :: send_count !< send count for MPI_GATHERV
	949
	950	INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_y !< length scale in y-direction
	951	INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z !< length scale in z-direction
	952	INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_y2 !< n_y*2
	953	INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z2 !< n_z*2
	954
	955	REAL(wp) :: nyz_inv !< inverse of number of grid points in yz-slice
	956	REAL(wp) :: rand_av !< average of random number
	957	REAL(wp) :: rand_sigma_inv !< inverse of stdev of random number
	958
	959	REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1) :: b_y !< filter func in y-dir
	960	REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1) :: b_z !< filter func in z-dir
	961	REAL(wp), DIMENSION(nzb_x:nzt_x+1,nysg:nyng) :: f_n_l !< local velocity seed
	962	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng) :: f_n !< velocity seed
	963	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rand_it !< random number
	964
	965
	966	!
	967	!-- Generate random numbers using a seed generated in stg_init.
	968	!-- The set of random numbers are modified to have an average of 0 and
	969	!-- unit variance.
	970	ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp:ny+mergp) )
	971
	972	rand_av = 0.0_wp
	973	rand_sigma_inv = 0.0_wp
	974	nyz_inv = 1.0_wp / REAL( ( nzt+1 - nzb+1 ) * ( ny+1 ), KIND=wp )
	975
	976	DO j = 0, ny
	977	DO k = nzb, nzt+1
	978	CALL RANDOM_NUMBER( rand_it(k,j) )
	979	rand_av = rand_av + rand_it(k,j)
	980	ENDDO
	981	ENDDO
	982
	983	rand_av = rand_av * nyz_inv
	984
	985	DO j = 0, ny
	986	DO k = nzb, nzt+1
	987	rand_it(k,j) = rand_it(k,j) - rand_av
	988	rand_sigma_inv = rand_sigma_inv + rand_it(k,j) ** 2
	989	ENDDO
	990	ENDDO
	991
	992	rand_sigma_inv = 1.0_wp / SQRT(rand_sigma_inv * nyz_inv)
	993
	994	DO j = 0, ny
	995	DO k = nzb, nzt+1
	996	rand_it(k,j) = rand_it(k,j) * rand_sigma_inv
	997	ENDDO
	998	ENDDO
	999
	1000	!
	1001	!-- Periodic fill of random number in space
	1002	DO j = 0, ny
	1003	DO k = 1, mergp
	1004	rand_it(nzb -k,j) = rand_it(nzt+2-k,j) ! bottom margin
	1005	rand_it(nzt+1+k,j) = rand_it(nzb+k-1,j) ! top margin
	1006	ENDDO
	1007	ENDDO
	1008	DO j = 1, mergp
	1009	DO k = nzb-mergp, nzt+1+mergp
	1010	rand_it(k, -j) = rand_it(k,ny-j+1) ! south margin
	1011	rand_it(k,ny+j) = rand_it(k, j-1) ! north margin
	1012	ENDDO
	1013	ENDDO
	1014
	1015	!
	1016	!-- Generate velocity seed following Eq.6 of Xie and Castro (2008)
	1017	n_y2 = n_y * 2
	1018	n_z2 = n_z * 2
	1019	f_n_l = 0.0_wp
	1020
	1021	DO j = nysg, nyng
	1022	DO k = nzb_x, nzt_x+1
	1023	DO jj = -n_y2(k), n_y2(k)
	1024	DO kk = -n_z2(k), n_z2(k)
	1025	f_n_l(k,j) = f_n_l(k,j) &
	1026	+ b_y(jj,k) * b_z(kk,k) * rand_it(k+kk,j+jj)
	1027	ENDDO
	1028	ENDDO
	1029	ENDDO
	1030	ENDDO
	1031
	1032	DEALLOCATE( rand_it )
	1033
	1034	!
	1035	!-- Gather velocity seeds of full subdomain
	1036	send_count = nzt_x - nzb_x + 1
	1037	IF ( nzt_x == nzt ) send_count = send_count + 1
	1038
	1039	#if defined( __parallel )
	1040	CALL MPI_GATHERV( f_n_l(nzb_x,nysg), send_count, stg_type_yz_small, &
	1041	f_n(nzb+1,nysg), recv_count, displs, stg_type_yz, &
	1042	id_inflow, comm1dx, ierr )
	1043	#else
	1044	f_n(nzb+1:nzt+1,nysg:nyng) = f_n_l(nzb_x:nzt_x+1,nysg:nyng)
	1045	#endif
	1046
	1047
	1048	END SUBROUTINE stg_generate_seed_yz
	1049
	1050	END MODULE synthetic_turbulence_generator_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |