source: palm/trunk/SOURCE/inflow_turbulence.f90 @ 1347

Last change on this file since 1347 was 1347, checked in by heinze, 7 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 7.1 KB
Line 
1 SUBROUTINE inflow_turbulence
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later 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 1997-2014 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: inflow_turbulence.f90 1347 2014-03-27 13:23:00Z heinze $
27!
28! 1346 2014-03-27 13:18:20Z heinze
29! Bugfix: REAL constants provided with KIND-attribute especially in call of
30! intrinsic function like MAX, MIN, SIGN
31!
32! 1320 2014-03-20 08:40:49Z raasch
33! ONLY-attribute added to USE-statements,
34! kind-parameters added to all INTEGER and REAL declaration statements,
35! kinds are defined in new module kinds,
36! revision history before 2012 removed,
37! comment fields (!:) to be used for variable explanations added to
38! all variable declaration statements
39!
40! 1092 2013-02-02 11:24:22Z raasch
41! unused variables removed
42!
43! 1036 2012-10-22 13:43:42Z raasch
44! code put under GPL (PALM 3.9)
45!
46! Initial version (2008/03/07)
47!
48! Description:
49! ------------
50! Imposing turbulence at the respective inflow using the turbulence
51! recycling method of Kataoka and Mizuno (2002).
52!------------------------------------------------------------------------------!
53
54    USE arrays_3d,                                                             &
55        ONLY:  e, inflow_damping_factor, mean_inflow_profiles, pt, u, v, w
56       
57    USE control_parameters,                                                    &
58        ONLY:  recycling_plane
59       
60    USE cpulog,                                                                &
61        ONLY:  cpu_log, log_point
62       
63    USE grid_variables,                                                        &
64        ONLY: 
65       
66    USE indices,                                                               &
67        ONLY:  nbgp, nxl, ny, nyn, nys, nyng, nysg, nzb, nzt
68       
69    USE kinds
70   
71    USE pegrid
72
73
74    IMPLICIT NONE
75
76    INTEGER(iwp) ::  i        !:
77    INTEGER(iwp) ::  j        !:
78    INTEGER(iwp) ::  k        !:
79    INTEGER(iwp) ::  l        !:
80    INTEGER(iwp) ::  ngp_ifd  !:
81    INTEGER(iwp) ::  ngp_pr   !:
82
83    REAL(wp), DIMENSION(nzb:nzt+1,5,nbgp)           ::                         &
84       avpr, avpr_l  !:
85    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,5,nbgp) ::                         &
86       inflow_dist   !:
87
88    CALL cpu_log( log_point(40), 'inflow_turbulence', 'start' )
89
90!
91!-- Carry out spanwise averaging in the recycling plane
92    avpr_l = 0.0
93    ngp_pr = ( nzt - nzb + 2 ) * 5 * nbgp
94    ngp_ifd = ngp_pr * ( nyn - nys + 1 + 2 * nbgp )
95
96!
97!-- First, local averaging within the recycling domain
98    i = recycling_plane
99
100#if defined( __parallel )
101    IF ( myidx == id_recycling )  THEN
102       
103       DO  l = 1, nbgp
104          DO  j = nys, nyn
105             DO  k = nzb, nzt + 1
106
107                avpr_l(k,1,l) = avpr_l(k,1,l) + u(k,j,i)
108                avpr_l(k,2,l) = avpr_l(k,2,l) + v(k,j,i)
109                avpr_l(k,3,l) = avpr_l(k,3,l) + w(k,j,i)
110                avpr_l(k,4,l) = avpr_l(k,4,l) + pt(k,j,i)
111                avpr_l(k,5,l) = avpr_l(k,5,l) + e(k,j,i)
112
113             ENDDO
114          ENDDO
115          i = i + 1
116       ENDDO
117
118    ENDIF
119!
120!-- Now, averaging over all PEs
121    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
122    CALL MPI_ALLREDUCE( avpr_l(nzb,1,1), avpr(nzb,1,1), ngp_pr, MPI_REAL, &
123                        MPI_SUM, comm2d, ierr )
124
125#else
126    DO  l = 1, nbgp
127       DO  j = nys, nyn
128          DO  k = nzb, nzt + 1
129
130             avpr_l(k,1,l) = avpr_l(k,1,l) + u(k,j,i)
131             avpr_l(k,2,l) = avpr_l(k,2,l) + v(k,j,i)
132             avpr_l(k,3,l) = avpr_l(k,3,l) + w(k,j,i)
133             avpr_l(k,4,l) = avpr_l(k,4,l) + pt(k,j,i)
134             avpr_l(k,5,l) = avpr_l(k,5,l) + e(k,j,i)
135
136          ENDDO
137       ENDDO
138       i = i + 1 
139    ENDDO
140   
141    avpr = avpr_l
142#endif
143
144    avpr = avpr / ( ny + 1 )
145!
146!-- Calculate the disturbances at the recycling plane
147    i = recycling_plane
148
149#if defined( __parallel )
150    IF ( myidx == id_recycling )  THEN
151       DO  l = 1, nbgp
152          DO  j = nysg, nyng
153             DO  k = nzb, nzt + 1
154
155                inflow_dist(k,j,1,l) = u(k,j,i+1) - avpr(k,1,l)
156                inflow_dist(k,j,2,l) = v(k,j,i)   - avpr(k,2,l)
157                inflow_dist(k,j,3,l) = w(k,j,i)   - avpr(k,3,l)
158                inflow_dist(k,j,4,l) = pt(k,j,i)  - avpr(k,4,l)
159                inflow_dist(k,j,5,l) = e(k,j,i)   - avpr(k,5,l)
160             
161            ENDDO
162          ENDDO
163          i = i + 1
164       ENDDO
165
166    ENDIF
167#else
168    DO  l = 1, nbgp
169       DO  j = nysg, nyng
170          DO  k = nzb, nzt+1
171
172             inflow_dist(k,j,1,l) = u(k,j,i+1) - avpr(k,1,l)
173             inflow_dist(k,j,2,l) = v(k,j,i)   - avpr(k,2,l)
174             inflow_dist(k,j,3,l) = w(k,j,i)   - avpr(k,3,l)
175             inflow_dist(k,j,4,l) = pt(k,j,i)  - avpr(k,4,l)
176             inflow_dist(k,j,5,l) = e(k,j,i)   - avpr(k,5,l)
177             
178          ENDDO
179       ENDDO
180       i = i + 1
181    ENDDO
182#endif
183
184!
185!-- For parallel runs, send the disturbances to the respective inflow PE
186#if defined( __parallel )
187    IF ( myidx == id_recycling  .AND.  myidx /= id_inflow )  THEN
188
189       CALL MPI_SEND( inflow_dist(nzb,nysg,1,1), ngp_ifd, MPI_REAL, &
190                      id_inflow, 1, comm1dx, ierr )
191
192    ELSEIF ( myidx /= id_recycling  .AND.  myidx == id_inflow )  THEN
193
194       inflow_dist = 0.0
195       CALL MPI_RECV( inflow_dist(nzb,nysg,1,1), ngp_ifd, MPI_REAL, &
196                      id_recycling, 1, comm1dx, status, ierr )
197
198    ENDIF
199#endif
200
201!
202!-- Add the disturbance at the inflow
203    IF ( nxl == 0 )  THEN
204
205       DO  j = nysg, nyng
206          DO  k = nzb, nzt + 1
207
208              u(k,j,-nbgp+1:0) = mean_inflow_profiles(k,1) + &
209                           inflow_dist(k,j,1,1:nbgp) * inflow_damping_factor(k)
210              v(k,j,-nbgp:-1)  = mean_inflow_profiles(k,2) + &
211                           inflow_dist(k,j,2,1:nbgp) * inflow_damping_factor(k)
212              w(k,j,-nbgp:-1)  =                             &
213                           inflow_dist(k,j,3,1:nbgp) * inflow_damping_factor(k)
214              pt(k,j,-nbgp:-1) = mean_inflow_profiles(k,4) + &
215                           inflow_dist(k,j,4,1:nbgp) * inflow_damping_factor(k)
216              e(k,j,-nbgp:-1)  = mean_inflow_profiles(k,5) + &
217                           inflow_dist(k,j,5,1:nbgp) * inflow_damping_factor(k)
218              e(k,j,-nbgp:-1)  = MAX( e(k,j,-nbgp:-1), 0.0_wp )
219
220          ENDDO
221       ENDDO
222
223    ENDIF
224
225    CALL cpu_log( log_point(40), 'inflow_turbulence', 'stop' )
226
227
228 END SUBROUTINE inflow_turbulence
Note: See TracBrowser for help on using the repository browser.