source: palm/trunk/SOURCE/lpm_init_sgs_tke.f90 @ 1360

Last change on this file since 1360 was 1360, checked in by hoffmann, 7 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 9.2 KB
Line 
1 SUBROUTINE lpm_init_sgs_tke
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: lpm_init_sgs_tke.f90 1360 2014-04-11 17:20:32Z hoffmann $
27!
28! 1359 2014-04-11 17:15:14Z hoffmann
29! New particle structure integrated.
30! Kind definition added to all floating point numbers.
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! comment fields (!:) to be used for variable explanations added to
37! all variable declaration statements
38!
39! 1036 2012-10-22 13:43:42Z raasch
40! code put under GPL (PALM 3.9)
41!
42! 849 2012-03-15 10:35:09Z raasch
43! initial revision (former part of advec_particles)
44!
45!
46! Description:
47! ------------
48! Calculates quantities required for considering the SGS velocity fluctuations
49! in the particle transport by a stochastic approach. The respective
50! quantities are: SGS-TKE gradients and horizontally averaged profiles of the
51! SGS TKE and the resolved-scale velocity variances.
52!------------------------------------------------------------------------------!
53
54    USE arrays_3d,                                                             &
55        ONLY:  de_dx, de_dy, de_dz, diss, e, u, v, w, zu
56
57    USE grid_variables,                                                        &
58        ONLY:  ddx, ddy
59
60    USE indices,                                                               &
61        ONLY:  nbgp, ngp_2dh_outer, nx, nxl, nxr, ny, nyn, nys, nz, nzb,       &
62                      nzb_s_inner, nzb_s_outer, nzt
63
64    USE kinds
65
66    USE particle_attributes,                                                   &
67        ONLY:  sgs_wfu_part, sgs_wfv_part, sgs_wfw_part
68
69    USE pegrid
70
71    USE statistics,                                                            &
72        ONLY:  flow_statistics_called, hom, sums, sums_l
73
74    IMPLICIT NONE
75
76    INTEGER(iwp) ::  i      !:
77    INTEGER(iwp) ::  j      !:
78    INTEGER(iwp) ::  k      !:
79
80
81!
82!-- TKE gradient along x and y
83    DO  i = nxl, nxr
84       DO  j = nys, nyn
85          DO  k = nzb, nzt+1
86
87             IF ( k <= nzb_s_inner(j,i-1)  .AND.  k > nzb_s_inner(j,i)  .AND.  &
88                  k  > nzb_s_inner(j,i+1) )                                    &
89             THEN
90                de_dx(k,j,i) = 2.0_wp * sgs_wfu_part *                         &
91                               ( e(k,j,i+1) - e(k,j,i) ) * ddx
92             ELSEIF ( k  > nzb_s_inner(j,i-1)  .AND.  k > nzb_s_inner(j,i)     &
93                      .AND.  k <= nzb_s_inner(j,i+1) )                         &
94             THEN
95                de_dx(k,j,i) = 2.0_wp * sgs_wfu_part *                         &
96                               ( e(k,j,i) - e(k,j,i-1) ) * ddx
97             ELSEIF ( k < nzb_s_inner(j,i)  .AND.  k < nzb_s_inner(j,i+1) )    &
98             THEN
99                de_dx(k,j,i) = 0.0_wp
100             ELSEIF ( k < nzb_s_inner(j,i-1)  .AND.  k < nzb_s_inner(j,i) )    &
101             THEN
102                de_dx(k,j,i) = 0.0_wp
103             ELSE
104                de_dx(k,j,i) = sgs_wfu_part * ( e(k,j,i+1) - e(k,j,i-1) ) * ddx
105             ENDIF
106
107             IF ( k <= nzb_s_inner(j-1,i)  .AND.  k > nzb_s_inner(j,i)  .AND.  &
108                  k  > nzb_s_inner(j+1,i) )                                    &
109             THEN
110                de_dy(k,j,i) = 2.0_wp * sgs_wfv_part *                         &
111                               ( e(k,j+1,i) - e(k,j,i) ) * ddy
112             ELSEIF ( k  > nzb_s_inner(j-1,i)  .AND.  k  > nzb_s_inner(j,i)    &
113                      .AND.  k <= nzb_s_inner(j+1,i) )                         &
114             THEN
115                de_dy(k,j,i) = 2.0_wp * sgs_wfv_part *                         &
116                               ( e(k,j,i) - e(k,j-1,i) ) * ddy
117             ELSEIF ( k < nzb_s_inner(j,i)  .AND.  k < nzb_s_inner(j+1,i) )    &
118             THEN
119                de_dy(k,j,i) = 0.0_wp
120             ELSEIF ( k < nzb_s_inner(j-1,i)  .AND.  k < nzb_s_inner(j,i) )    &
121             THEN
122                de_dy(k,j,i) = 0.0_wp
123             ELSE
124                de_dy(k,j,i) = sgs_wfv_part * ( e(k,j+1,i) - e(k,j-1,i) ) * ddy
125             ENDIF
126
127          ENDDO
128       ENDDO
129    ENDDO
130
131!
132!-- TKE gradient along z, including bottom and top boundary conditions
133    DO  i = nxl, nxr
134       DO  j = nys, nyn
135
136          DO  k = nzb_s_inner(j,i)+2, nzt-1
137             de_dz(k,j,i)  = 2.0_wp * sgs_wfw_part *                           &
138                             ( e(k+1,j,i) - e(k-1,j,i) ) / ( zu(k+1)-zu(k-1) )
139          ENDDO
140
141          k = nzb_s_inner(j,i)
142          de_dz(nzb:k,j,i) = 0.0_wp
143          de_dz(k+1,j,i)   = 2.0_wp * sgs_wfw_part *                           &
144                             ( e(k+2,j,i) - e(k+1,j,i) ) / ( zu(k+2) - zu(k+1) )
145          de_dz(nzt,j,i)   = 0.0_wp
146          de_dz(nzt+1,j,i) = 0.0_wp
147       ENDDO
148    ENDDO
149
150
151!
152!-- Lateral boundary conditions
153    CALL exchange_horiz( de_dx, nbgp )
154    CALL exchange_horiz( de_dy, nbgp )
155    CALL exchange_horiz( de_dz, nbgp )
156    CALL exchange_horiz( diss, nbgp  )
157
158
159!
160!-- Calculate the horizontally averaged profiles of SGS TKE and resolved
161!-- velocity variances (they may have been already calculated in routine
162!-- flow_statistics).
163    IF ( .NOT. flow_statistics_called )  THEN
164
165!
166!--    First calculate horizontally averaged profiles of the horizontal
167!--    velocities.
168       sums_l(:,1,0) = 0.0_wp
169       sums_l(:,2,0) = 0.0_wp
170
171       DO  i = nxl, nxr
172          DO  j =  nys, nyn
173             DO  k = nzb_s_outer(j,i), nzt+1
174                sums_l(k,1,0)  = sums_l(k,1,0)  + u(k,j,i)
175                sums_l(k,2,0)  = sums_l(k,2,0)  + v(k,j,i)
176             ENDDO
177          ENDDO
178       ENDDO
179
180#if defined( __parallel )
181!
182!--    Compute total sum from local sums
183       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
184       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, &
185                           MPI_REAL, MPI_SUM, comm2d, ierr )
186       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
187       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, &
188                              MPI_REAL, MPI_SUM, comm2d, ierr )
189#else
190       sums(:,1) = sums_l(:,1,0)
191       sums(:,2) = sums_l(:,2,0)
192#endif
193
194!
195!--    Final values are obtained by division by the total number of grid
196!--    points used for the summation.
197       hom(:,1,1,0) = sums(:,1) / ngp_2dh_outer(:,0)   ! u
198       hom(:,1,2,0) = sums(:,2) / ngp_2dh_outer(:,0)   ! v
199
200!
201!--    Now calculate the profiles of SGS TKE and the resolved-scale
202!--    velocity variances
203       sums_l(:,8,0)  = 0.0_wp
204       sums_l(:,30,0) = 0.0_wp
205       sums_l(:,31,0) = 0.0_wp
206       sums_l(:,32,0) = 0.0_wp
207       DO  i = nxl, nxr
208          DO  j = nys, nyn
209             DO  k = nzb_s_outer(j,i), nzt+1
210                sums_l(k,8,0)  = sums_l(k,8,0)  + e(k,j,i)
211                sums_l(k,30,0) = sums_l(k,30,0) + ( u(k,j,i) - hom(k,1,1,0) )**2
212                sums_l(k,31,0) = sums_l(k,31,0) + ( v(k,j,i) - hom(k,1,2,0) )**2
213                sums_l(k,32,0) = sums_l(k,32,0) + w(k,j,i)**2
214             ENDDO
215          ENDDO
216       ENDDO
217
218#if defined( __parallel )
219!
220!--    Compute total sum from local sums
221       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
222       CALL MPI_ALLREDUCE( sums_l(nzb,8,0), sums(nzb,8), nzt+2-nzb, &
223                           MPI_REAL, MPI_SUM, comm2d, ierr )
224       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
225       CALL MPI_ALLREDUCE( sums_l(nzb,30,0), sums(nzb,30), nzt+2-nzb, &
226                           MPI_REAL, MPI_SUM, comm2d, ierr )
227       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
228       CALL MPI_ALLREDUCE( sums_l(nzb,31,0), sums(nzb,31), nzt+2-nzb, &
229                           MPI_REAL, MPI_SUM, comm2d, ierr )
230       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
231       CALL MPI_ALLREDUCE( sums_l(nzb,32,0), sums(nzb,32), nzt+2-nzb, &
232                           MPI_REAL, MPI_SUM, comm2d, ierr )
233
234#else
235       sums(:,8)  = sums_l(:,8,0)
236       sums(:,30) = sums_l(:,30,0)
237       sums(:,31) = sums_l(:,31,0)
238       sums(:,32) = sums_l(:,32,0)
239#endif
240
241!
242!--    Final values are obtained by division by the total number of grid
243!--    points used for the summation.
244       hom(:,1,8,0)  = sums(:,8)  / ngp_2dh_outer(:,0)   ! e
245       hom(:,1,30,0) = sums(:,30) / ngp_2dh_outer(:,0)   ! u*2
246       hom(:,1,31,0) = sums(:,31) / ngp_2dh_outer(:,0)   ! v*2
247       hom(:,1,32,0) = sums(:,32) / ngp_2dh_outer(:,0)   ! w*2
248
249    ENDIF
250
251 END SUBROUTINE lpm_init_sgs_tke
Note: See TracBrowser for help on using the repository browser.