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

Last change on this file since 1930 was 1930, checked in by suehring, 5 years ago

last commit documented

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