Home

Context Navigation

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

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