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

Last change on this file since 2036 was 2001, checked in by knoop, 8 years ago

last commit documented

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