Home

Context Navigation

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

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