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

Last change on this file since 2332 was 2259, checked in by gronemeier, 7 years ago

Implemented synthetic turbulence generator

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