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

Last change on this file since 668 was 668, checked in by suehring, 11 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 6.3 KB
Line 
1 SUBROUTINE inflow_turbulence
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6!
7! Former revisions:
8! -----------------
9! $Id: inflow_turbulence.f90 668 2010-12-23 13:22:58Z suehring $
10!
11! 667 2010-12-23 12:06:00Z suehring/gryschka
12! Using nbgp recycling planes for a better resolution of the turbulent flow
13! near the inflow.
14!
15! 622 2010-12-10 08:08:13Z raasch
16! optional barriers included in order to speed up collective operations
17!
18! 222 2009-01-12 16:04:16Z letzel
19! Bugfix for nonparallel execution
20!
21! Initial version (2008/03/07)
22!
23! Description:
24! ------------
25! Imposing turbulence at the respective inflow using the turbulence
26! recycling method of Kataoka and Mizuno (2002).
27!------------------------------------------------------------------------------!
28
29    USE arrays_3d
30    USE control_parameters
31    USE cpulog
32    USE grid_variables
33    USE indices
34    USE interfaces
35    USE pegrid
36
37
38    IMPLICIT NONE
39
40    INTEGER ::  i, imax, j, k, l, ngp_ifd, ngp_pr
41
42    REAL, DIMENSION(1:2) ::  volume_flow_l, volume_flow_offset
43    REAL, DIMENSION(nzb:nzt+1,5,nbgp) ::  avpr, avpr_l
44    REAL, DIMENSION(nzb:nzt+1,nysg:nyng,5,nbgp) ::  inflow_dist
45
46    CALL cpu_log( log_point(40), 'inflow_turbulence', 'start' )
47
48!
49!-- Carry out spanwise averaging in the recycling plane
50    avpr_l = 0.0
51    ngp_pr = ( nzt - nzb + 2 ) * 5 * nbgp
52    ngp_ifd = ngp_pr * ( nyn - nys + 1 + 2 * nbgp )
53
54!
55!-- First, local averaging within the recycling domain
56
57    i = recycling_plane
58
59#if defined( __parallel )
60    IF ( myidx == id_recycling )  THEN
61       
62       DO  l = 1, nbgp
63          DO  j = nys, nyn
64             DO  k = nzb, nzt + 1
65
66                avpr_l(k,1,l) = avpr_l(k,1,l) + u(k,j,i)
67                avpr_l(k,2,l) = avpr_l(k,2,l) + v(k,j,i)
68                avpr_l(k,3,l) = avpr_l(k,3,l) + w(k,j,i)
69                avpr_l(k,4,l) = avpr_l(k,4,l) + pt(k,j,i)
70                avpr_l(k,5,l) = avpr_l(k,5,l) + e(k,j,i)
71
72             ENDDO
73          ENDDO
74          i = i + 1
75       ENDDO
76
77    ENDIF
78!
79!-- Now, averaging over all PEs
80    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
81    CALL MPI_ALLREDUCE( avpr_l(nzb,1,1), avpr(nzb,1,1), ngp_pr,  & 
82                    MPI_REAL, MPI_SUM, comm2d, ierr )
83
84#else
85    DO  l = 1, nbgp
86       DO  j = nys, nyn
87          DO  k = nzb, nzt + 1
88
89             avpr_l(k,1,l) = avpr_l(k,1,l) + u(k,j,i)
90             avpr_l(k,2,l) = avpr_l(k,2,l) + v(k,j,i)
91             avpr_l(k,3,l) = avpr_l(k,3,l) + w(k,j,i)
92             avpr_l(k,4,l) = avpr_l(k,4,l) + pt(k,j,i)
93             avpr_l(k,5,l) = avpr_l(k,5,l) + e(k,j,i)
94
95          ENDDO
96       ENDDO
97       i = i + 1 
98    ENDDO
99   
100    avpr = avpr_l
101#endif
102
103    avpr = avpr / ( ny + 1 )
104!
105!-- Calculate the disturbances at the recycling plane
106    i = recycling_plane
107
108#if defined( __parallel )
109    IF ( myidx == id_recycling )  THEN
110       DO  l = 1, nbgp
111          DO  j = nysg, nyng
112             DO  k = nzb, nzt + 1
113
114                inflow_dist(k,j,1,l) = u(k,j,i+1) - avpr(k,1,l)
115                inflow_dist(k,j,2,l) = v(k,j,i)   - avpr(k,2,l)
116                inflow_dist(k,j,3,l) = w(k,j,i)   - avpr(k,3,l)
117                inflow_dist(k,j,4,l) = pt(k,j,i)  - avpr(k,4,l)
118                inflow_dist(k,j,5,l) = e(k,j,i)   - avpr(k,5,l)
119             
120            ENDDO
121          ENDDO
122          i = i + 1
123       ENDDO
124
125    ENDIF
126#else
127    DO  l = 1, nbgp
128       DO  j = nysg, nyng
129          DO  k = nzb, nzt+1
130
131             inflow_dist(k,j,1,l) = u(k,j,i+1) - avpr(k,1,l)
132             inflow_dist(k,j,2,l) = v(k,j,i)   - avpr(k,2,l)
133             inflow_dist(k,j,3,l) = w(k,j,i)   - avpr(k,3,l)
134             inflow_dist(k,j,4,l) = pt(k,j,i)  - avpr(k,4,l)
135             inflow_dist(k,j,5,l) = e(k,j,i)   - avpr(k,5,l)
136             
137          ENDDO
138       ENDDO
139       i = i + 1
140    ENDDO
141#endif
142
143!
144!-- For parallel runs, send the disturbances to the respective inflow PE
145#if defined( __parallel )
146    IF ( myidx == id_recycling  .AND.  myidx /= id_inflow )  THEN
147
148       CALL MPI_SEND( inflow_dist(nzb,nysg,1,1), ngp_ifd, MPI_REAL, &
149                      id_inflow, 1, comm1dx, ierr )
150
151    ELSEIF ( myidx /= id_recycling  .AND.  myidx == id_inflow )  THEN
152
153       inflow_dist = 0.0
154       CALL MPI_RECV( inflow_dist(nzb,nysg,1,1), ngp_ifd, MPI_REAL, &
155                      id_recycling, 1, comm1dx, status, ierr )
156
157    ENDIF
158#endif
159
160!
161!-- Add the disturbance at the inflow
162    IF ( nxl == 0 )  THEN
163
164       DO  j = nysg, nyng
165          DO  k = nzb, nzt + 1
166
167              u(k,j,-nbgp+1:0)   = mean_inflow_profiles(k,1) + &
168                           inflow_dist(k,j,1,1:nbgp) * inflow_damping_factor(k)
169              v(k,j,-nbgp:-1)  = mean_inflow_profiles(k,2) + &
170                           inflow_dist(k,j,2,1:nbgp) * inflow_damping_factor(k)
171              w(k,j,-nbgp:-1)  = inflow_dist(k,j,3,1:nbgp) * inflow_damping_factor(k)
172              pt(k,j,-nbgp:-1) = mean_inflow_profiles(k,4) + &
173                           inflow_dist(k,j,4,1:nbgp) * inflow_damping_factor(k)
174              e(k,j,-nbgp:-1)  = mean_inflow_profiles(k,5) + &
175                           inflow_dist(k,j,5,1:nbgp) * inflow_damping_factor(k)
176              e(k,j,-nbgp:-1)  = MAX( e(k,j,-nbgp:-1), 0.0 )
177
178          ENDDO
179       ENDDO
180
181    ENDIF
182
183!
184!-- Conserve the volume flow at the inflow in order to avoid generation of
185!-- waves in the stable layer
186!    IF ( conserve_volume_flow  .AND.  inflow_l )  THEN
187
188!       volume_flow(1)   = 0.0
189!       volume_flow_l(1) = 0.0
190
191!       i = 0
192
193!       DO  j = nys, nyn
194!
195!--       Sum up the volume flow through the south/north boundary
196!          DO  k = nzb_2d(j,i) + 1, nzt
197!             volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzu(k)
198!          ENDDO
199!       ENDDO
200
201#if defined( __parallel )   
202!       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
203!       CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, &
204!                           MPI_SUM, comm1dy, ierr )   
205#else
206!       volume_flow = volume_flow_l 
207#endif
208!       volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) )    &
209!                               / volume_flow_area(1)
210
211!       DO  j = nys-1, nyn+1
212!          DO  k = nzb_v_inner(j,i) + 1, nzt
213!             u(k,j,i) = u(k,j,i) + volume_flow_offset(1)
214!          ENDDO
215!       ENDDO
216
217!    ENDIF
218
219    CALL cpu_log( log_point(40), 'inflow_turbulence', 'stop' )
220
221
222 END SUBROUTINE inflow_turbulence
Note: See TracBrowser for help on using the repository browser.