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

Last change on this file since 2638 was 2576, checked in by Giersch, 6 years ago

Bugfixes for restart runs

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