Home

Context Navigation

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

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