Home

Context Navigation

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

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