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

Last change on this file since 1818 was 1818, checked in by maronga, 5 years ago

last commit documented / copyright update

  • Property svn:keywords set to Id
File size: 9.3 KB
Line 
1!> @file lpm_init_sgs_tke.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2016 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! ------------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: lpm_init_sgs_tke.f90 1818 2016-04-06 15:53:27Z maronga $
26!
27! 1682 2015-10-07 23:56:08Z knoop
28! Code annotations made doxygen readable
29!
30! 1359 2014-04-11 17:15:14Z hoffmann
31! New particle structure integrated.
32! Kind definition added to all floating point numbers.
33!
34! 1320 2014-03-20 08:40:49Z raasch
35! ONLY-attribute added to USE-statements,
36! kind-parameters added to all INTEGER and REAL declaration statements,
37! kinds are defined in new module kinds,
38! comment fields (!:) to be used for variable explanations added to
39! all variable declaration statements
40!
41! 1036 2012-10-22 13:43:42Z raasch
42! code put under GPL (PALM 3.9)
43!
44! 849 2012-03-15 10:35:09Z raasch
45! initial revision (former part of advec_particles)
46!
47!
48! Description:
49! ------------
50!> Calculates quantities required for considering the SGS velocity fluctuations
51!> in the particle transport by a stochastic approach. The respective
52!> quantities are: SGS-TKE gradients and horizontally averaged profiles of the
53!> SGS TKE and the resolved-scale velocity variances.
54!------------------------------------------------------------------------------!
55 SUBROUTINE lpm_init_sgs_tke
56 
57
58    USE arrays_3d,                                                             &
59        ONLY:  de_dx, de_dy, de_dz, diss, e, u, v, w, zu
60
61    USE grid_variables,                                                        &
62        ONLY:  ddx, ddy
63
64    USE indices,                                                               &
65        ONLY:  nbgp, ngp_2dh_outer, nx, nxl, nxr, ny, nyn, nys, nz, nzb,       &
66                      nzb_s_inner, nzb_s_outer, nzt
67
68    USE kinds
69
70    USE particle_attributes,                                                   &
71        ONLY:  sgs_wfu_part, sgs_wfv_part, sgs_wfw_part
72
73    USE pegrid
74
75    USE statistics,                                                            &
76        ONLY:  flow_statistics_called, hom, sums, sums_l
77
78    IMPLICIT NONE
79
80    INTEGER(iwp) ::  i      !<
81    INTEGER(iwp) ::  j      !<
82    INTEGER(iwp) ::  k      !<
83
84
85!
86!-- TKE gradient along x and y
87    DO  i = nxl, nxr
88       DO  j = nys, nyn
89          DO  k = nzb, nzt+1
90
91             IF ( k <= nzb_s_inner(j,i-1)  .AND.  k > nzb_s_inner(j,i)  .AND.  &
92                  k  > nzb_s_inner(j,i+1) )                                    &
93             THEN
94                de_dx(k,j,i) = 2.0_wp * sgs_wfu_part *                         &
95                               ( e(k,j,i+1) - e(k,j,i) ) * ddx
96             ELSEIF ( k  > nzb_s_inner(j,i-1)  .AND.  k > nzb_s_inner(j,i)     &
97                      .AND.  k <= nzb_s_inner(j,i+1) )                         &
98             THEN
99                de_dx(k,j,i) = 2.0_wp * sgs_wfu_part *                         &
100                               ( e(k,j,i) - e(k,j,i-1) ) * ddx
101             ELSEIF ( k < nzb_s_inner(j,i)  .AND.  k < nzb_s_inner(j,i+1) )    &
102             THEN
103                de_dx(k,j,i) = 0.0_wp
104             ELSEIF ( k < nzb_s_inner(j,i-1)  .AND.  k < nzb_s_inner(j,i) )    &
105             THEN
106                de_dx(k,j,i) = 0.0_wp
107             ELSE
108                de_dx(k,j,i) = sgs_wfu_part * ( e(k,j,i+1) - e(k,j,i-1) ) * ddx
109             ENDIF
110
111             IF ( k <= nzb_s_inner(j-1,i)  .AND.  k > nzb_s_inner(j,i)  .AND.  &
112                  k  > nzb_s_inner(j+1,i) )                                    &
113             THEN
114                de_dy(k,j,i) = 2.0_wp * sgs_wfv_part *                         &
115                               ( e(k,j+1,i) - e(k,j,i) ) * ddy
116             ELSEIF ( k  > nzb_s_inner(j-1,i)  .AND.  k  > nzb_s_inner(j,i)    &
117                      .AND.  k <= nzb_s_inner(j+1,i) )                         &
118             THEN
119                de_dy(k,j,i) = 2.0_wp * sgs_wfv_part *                         &
120                               ( e(k,j,i) - e(k,j-1,i) ) * ddy
121             ELSEIF ( k < nzb_s_inner(j,i)  .AND.  k < nzb_s_inner(j+1,i) )    &
122             THEN
123                de_dy(k,j,i) = 0.0_wp
124             ELSEIF ( k < nzb_s_inner(j-1,i)  .AND.  k < nzb_s_inner(j,i) )    &
125             THEN
126                de_dy(k,j,i) = 0.0_wp
127             ELSE
128                de_dy(k,j,i) = sgs_wfv_part * ( e(k,j+1,i) - e(k,j-1,i) ) * ddy
129             ENDIF
130
131          ENDDO
132       ENDDO
133    ENDDO
134
135!
136!-- TKE gradient along z, including bottom and top boundary conditions
137    DO  i = nxl, nxr
138       DO  j = nys, nyn
139
140          DO  k = nzb_s_inner(j,i)+2, nzt-1
141             de_dz(k,j,i)  = 2.0_wp * sgs_wfw_part *                           &
142                             ( e(k+1,j,i) - e(k-1,j,i) ) / ( zu(k+1)-zu(k-1) )
143          ENDDO
144
145          k = nzb_s_inner(j,i)
146          de_dz(nzb:k,j,i) = 0.0_wp
147          de_dz(k+1,j,i)   = 2.0_wp * sgs_wfw_part *                           &
148                             ( e(k+2,j,i) - e(k+1,j,i) ) / ( zu(k+2) - zu(k+1) )
149          de_dz(nzt,j,i)   = 0.0_wp
150          de_dz(nzt+1,j,i) = 0.0_wp
151       ENDDO
152    ENDDO
153
154
155!
156!-- Lateral boundary conditions
157    CALL exchange_horiz( de_dx, nbgp )
158    CALL exchange_horiz( de_dy, nbgp )
159    CALL exchange_horiz( de_dz, nbgp )
160    CALL exchange_horiz( diss, nbgp  )
161
162
163!
164!-- Calculate the horizontally averaged profiles of SGS TKE and resolved
165!-- velocity variances (they may have been already calculated in routine
166!-- flow_statistics).
167    IF ( .NOT. flow_statistics_called )  THEN
168
169!
170!--    First calculate horizontally averaged profiles of the horizontal
171!--    velocities.
172       sums_l(:,1,0) = 0.0_wp
173       sums_l(:,2,0) = 0.0_wp
174
175       DO  i = nxl, nxr
176          DO  j =  nys, nyn
177             DO  k = nzb_s_outer(j,i), nzt+1
178                sums_l(k,1,0)  = sums_l(k,1,0)  + u(k,j,i)
179                sums_l(k,2,0)  = sums_l(k,2,0)  + v(k,j,i)
180             ENDDO
181          ENDDO
182       ENDDO
183
184#if defined( __parallel )
185!
186!--    Compute total sum from local sums
187       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
188       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, &
189                           MPI_REAL, MPI_SUM, comm2d, ierr )
190       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
191       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, &
192                              MPI_REAL, MPI_SUM, comm2d, ierr )
193#else
194       sums(:,1) = sums_l(:,1,0)
195       sums(:,2) = sums_l(:,2,0)
196#endif
197
198!
199!--    Final values are obtained by division by the total number of grid
200!--    points used for the summation.
201       hom(:,1,1,0) = sums(:,1) / ngp_2dh_outer(:,0)   ! u
202       hom(:,1,2,0) = sums(:,2) / ngp_2dh_outer(:,0)   ! v
203
204!
205!--    Now calculate the profiles of SGS TKE and the resolved-scale
206!--    velocity variances
207       sums_l(:,8,0)  = 0.0_wp
208       sums_l(:,30,0) = 0.0_wp
209       sums_l(:,31,0) = 0.0_wp
210       sums_l(:,32,0) = 0.0_wp
211       DO  i = nxl, nxr
212          DO  j = nys, nyn
213             DO  k = nzb_s_outer(j,i), nzt+1
214                sums_l(k,8,0)  = sums_l(k,8,0)  + e(k,j,i)
215                sums_l(k,30,0) = sums_l(k,30,0) + ( u(k,j,i) - hom(k,1,1,0) )**2
216                sums_l(k,31,0) = sums_l(k,31,0) + ( v(k,j,i) - hom(k,1,2,0) )**2
217                sums_l(k,32,0) = sums_l(k,32,0) + w(k,j,i)**2
218             ENDDO
219          ENDDO
220       ENDDO
221
222#if defined( __parallel )
223!
224!--    Compute total sum from local sums
225       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
226       CALL MPI_ALLREDUCE( sums_l(nzb,8,0), sums(nzb,8), nzt+2-nzb, &
227                           MPI_REAL, MPI_SUM, comm2d, ierr )
228       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
229       CALL MPI_ALLREDUCE( sums_l(nzb,30,0), sums(nzb,30), nzt+2-nzb, &
230                           MPI_REAL, MPI_SUM, comm2d, ierr )
231       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
232       CALL MPI_ALLREDUCE( sums_l(nzb,31,0), sums(nzb,31), nzt+2-nzb, &
233                           MPI_REAL, MPI_SUM, comm2d, ierr )
234       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
235       CALL MPI_ALLREDUCE( sums_l(nzb,32,0), sums(nzb,32), nzt+2-nzb, &
236                           MPI_REAL, MPI_SUM, comm2d, ierr )
237
238#else
239       sums(:,8)  = sums_l(:,8,0)
240       sums(:,30) = sums_l(:,30,0)
241       sums(:,31) = sums_l(:,31,0)
242       sums(:,32) = sums_l(:,32,0)
243#endif
244
245!
246!--    Final values are obtained by division by the total number of grid
247!--    points used for the summation.
248       hom(:,1,8,0)  = sums(:,8)  / ngp_2dh_outer(:,0)   ! e
249       hom(:,1,30,0) = sums(:,30) / ngp_2dh_outer(:,0)   ! u*2
250       hom(:,1,31,0) = sums(:,31) / ngp_2dh_outer(:,0)   ! v*2
251       hom(:,1,32,0) = sums(:,32) / ngp_2dh_outer(:,0)   ! w*2
252
253    ENDIF
254
255 END SUBROUTINE lpm_init_sgs_tke
Note: See TracBrowser for help on using the repository browser.