Home

Context Navigation

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

Last change on this file since 2566 was 2563, checked in by Giersch, 7 years ago
Restart runs with the usage of the wind turbine model are possible now. Further small at reading/writing restart data
Property svn:keywords set to `Id`
File size: 34.2 KB

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