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

Last change on this file since 2930 was 2894, checked in by Giersch, 6 years ago

Reading/Writing? data in case of restart runs revised

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