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

Last change on this file since 2232 was 2232, checked in by suehring, 4 years ago

Adjustments according new topography and surface-modelling concept implemented

  • Property svn:keywords set to Id
File size: 11.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
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-2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22! Adjustments according to new topography realization
23!
24! Former revisions:
25! -----------------
26! $Id: lpm_init_sgs_tke.f90 2232 2017-05-30 17:47:52Z 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,                   &
76               nzb_s_outer, nzt, wall_flags_0
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    USE surface_mod,                                                           &
89        ONLY:  bc_h
90
91    IMPLICIT NONE
92
93    INTEGER(iwp) ::  i      !< index variable along x
94    INTEGER(iwp) ::  j      !< index variable along y
95    INTEGER(iwp) ::  k      !< index variable along z
96    INTEGER(iwp) ::  m      !< running index for the surface elements
97
98    REAL(wp) ::  flag1      !< flag to mask topography
99
100!
101!-- TKE gradient along x and y
102    DO  i = nxl, nxr
103       DO  j = nys, nyn
104          DO  k = nzb, nzt+1
105
106             IF ( .NOT. BTEST( wall_flags_0(k,j,i-1), 0 )  .AND.               &
107                        BTEST( wall_flags_0(k,j,i), 0   )  .AND.               &
108                        BTEST( wall_flags_0(k,j,i+1), 0 ) )                    &
109             THEN
110                de_dx(k,j,i) = 2.0_wp * sgs_wf_part *                          &
111                               ( e(k,j,i+1) - e(k,j,i) ) * ddx
112             ELSEIF ( BTEST( wall_flags_0(k,j,i-1), 0 )  .AND.                 &
113                      BTEST( wall_flags_0(k,j,i), 0   )  .AND.                 &
114                .NOT. BTEST( wall_flags_0(k,j,i+1), 0 ) )                      &
115             THEN
116                de_dx(k,j,i) = 2.0_wp * sgs_wf_part *                          &
117                               ( e(k,j,i) - e(k,j,i-1) ) * ddx
118             ELSEIF ( .NOT. BTEST( wall_flags_0(k,j,i), 22   )  .AND.          &
119                      .NOT. BTEST( wall_flags_0(k,j,i+1), 22 ) )               &   
120             THEN
121                de_dx(k,j,i) = 0.0_wp
122             ELSEIF ( .NOT. BTEST( wall_flags_0(k,j,i-1), 22 )  .AND.          &
123                      .NOT. BTEST( wall_flags_0(k,j,i), 22   ) )               &
124             THEN
125                de_dx(k,j,i) = 0.0_wp
126             ELSE
127                de_dx(k,j,i) = sgs_wf_part * ( e(k,j,i+1) - e(k,j,i-1) ) * ddx
128             ENDIF
129
130             IF ( .NOT. BTEST( wall_flags_0(k,j-1,i), 0 )  .AND.               &
131                        BTEST( wall_flags_0(k,j,i), 0   )  .AND.               &
132                        BTEST( wall_flags_0(k,j+1,i), 0 ) )                    &
133             THEN
134                de_dy(k,j,i) = 2.0_wp * sgs_wf_part *                          &
135                               ( e(k,j+1,i) - e(k,j,i) ) * ddy
136             ELSEIF ( BTEST( wall_flags_0(k,j-1,i), 0 )  .AND.                 &
137                      BTEST( wall_flags_0(k,j,i), 0   )  .AND.                 &
138                .NOT. BTEST( wall_flags_0(k,j+1,i), 0 ) )                      &
139             THEN
140                de_dy(k,j,i) = 2.0_wp * sgs_wf_part *                          &
141                               ( e(k,j,i) - e(k,j-1,i) ) * ddy
142             ELSEIF ( .NOT. BTEST( wall_flags_0(k,j,i), 22   )  .AND.          &
143                      .NOT. BTEST( wall_flags_0(k,j+1,i), 22 ) )               &   
144             THEN
145                de_dy(k,j,i) = 0.0_wp
146             ELSEIF ( .NOT. BTEST( wall_flags_0(k,j-1,i), 22 )  .AND.          &
147                      .NOT. BTEST( wall_flags_0(k,j,i), 22   ) )               &
148             THEN
149                de_dy(k,j,i) = 0.0_wp
150             ELSE
151                de_dy(k,j,i) = sgs_wf_part * ( e(k,j+1,i) - e(k,j-1,i) ) * ddy
152             ENDIF
153
154          ENDDO
155       ENDDO
156    ENDDO
157
158!
159!-- TKE gradient along z at topograhy and  including bottom and top boundary conditions
160    DO  i = nxl, nxr
161       DO  j = nys, nyn
162          DO  k = nzb+1, nzt-1
163!
164!--          Flag to mask topography
165             flag1 = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0  ) )
166
167             de_dz(k,j,i) = 2.0_wp * sgs_wf_part *                             &
168                           ( e(k+1,j,i) - e(k-1,j,i) ) / ( zu(k+1) - zu(k-1) ) &
169                                                 * flag1
170          ENDDO
171!
172!--       upward-facing surfaces
173          DO  m = bc_h(0)%start_index(j,i), bc_h(0)%end_index(j,i)
174             k            = bc_h(0)%k(m)
175             de_dz(k,j,i) = 2.0_wp * sgs_wf_part *                             &
176                           ( e(k+1,j,i) - e(k,j,i)   ) / ( zu(k+1) - zu(k) )
177          ENDDO
178!
179!--       downward-facing surfaces
180          DO  m = bc_h(1)%start_index(j,i), bc_h(1)%end_index(j,i)
181             k            = bc_h(1)%k(m)
182             de_dz(k,j,i) = 2.0_wp * sgs_wf_part *                             &
183                           ( e(k,j,i) - e(k-1,j,i)   ) / ( zu(k) - zu(k-1) )
184          ENDDO
185
186          de_dz(nzb,j,i)   = 0.0_wp
187          de_dz(nzt,j,i)   = 0.0_wp
188          de_dz(nzt+1,j,i) = 0.0_wp
189       ENDDO
190    ENDDO
191!
192!-- Ghost point exchange
193    CALL exchange_horiz( de_dx, nbgp )
194    CALL exchange_horiz( de_dy, nbgp )
195    CALL exchange_horiz( de_dz, nbgp )
196    CALL exchange_horiz( diss, nbgp  )
197
198
199!
200!-- Calculate the horizontally averaged profiles of SGS TKE and resolved
201!-- velocity variances (they may have been already calculated in routine
202!-- flow_statistics).
203    IF ( .NOT. flow_statistics_called )  THEN
204
205!
206!--    First calculate horizontally averaged profiles of the horizontal
207!--    velocities.
208       sums_l(:,1,0) = 0.0_wp
209       sums_l(:,2,0) = 0.0_wp
210
211       DO  i = nxl, nxr
212          DO  j =  nys, nyn
213             DO  k = nzb, nzt+1
214!
215!--             Flag indicate nzb_s_outer
216                flag1 = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 24 ) )
217
218                sums_l(k,1,0)  = sums_l(k,1,0)  + u(k,j,i) * flag1
219                sums_l(k,2,0)  = sums_l(k,2,0)  + v(k,j,i) * flag1
220             ENDDO
221          ENDDO
222       ENDDO
223
224#if defined( __parallel )
225!
226!--    Compute total sum from local sums
227       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
228       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), 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,2,0), sums(nzb,2), nzt+2-nzb, &
232                              MPI_REAL, MPI_SUM, comm2d, ierr )
233#else
234       sums(:,1) = sums_l(:,1,0)
235       sums(:,2) = sums_l(:,2,0)
236#endif
237
238!
239!--    Final values are obtained by division by the total number of grid
240!--    points used for the summation.
241       hom(:,1,1,0) = sums(:,1) / ngp_2dh_outer(:,0)   ! u
242       hom(:,1,2,0) = sums(:,2) / ngp_2dh_outer(:,0)   ! v
243
244!
245!--    Now calculate the profiles of SGS TKE and the resolved-scale
246!--    velocity variances
247       sums_l(:,8,0)  = 0.0_wp
248       sums_l(:,30,0) = 0.0_wp
249       sums_l(:,31,0) = 0.0_wp
250       sums_l(:,32,0) = 0.0_wp
251       DO  i = nxl, nxr
252          DO  j = nys, nyn
253             DO  k = nzb, nzt+1
254!
255!--             Flag indicate nzb_s_outer
256                flag1 = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 24 ) )
257
258                sums_l(k,8,0)  = sums_l(k,8,0)  + e(k,j,i)                       * flag1
259                sums_l(k,30,0) = sums_l(k,30,0) + ( u(k,j,i) - hom(k,1,1,0) )**2 * flag1
260                sums_l(k,31,0) = sums_l(k,31,0) + ( v(k,j,i) - hom(k,1,2,0) )**2 * flag1
261                sums_l(k,32,0) = sums_l(k,32,0) + w(k,j,i)**2                    * flag1
262             ENDDO
263          ENDDO
264       ENDDO
265
266#if defined( __parallel )
267!
268!--    Compute total sum from local sums
269       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
270       CALL MPI_ALLREDUCE( sums_l(nzb,8,0), sums(nzb,8), nzt+2-nzb, &
271                           MPI_REAL, MPI_SUM, comm2d, ierr )
272       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
273       CALL MPI_ALLREDUCE( sums_l(nzb,30,0), sums(nzb,30), nzt+2-nzb, &
274                           MPI_REAL, MPI_SUM, comm2d, ierr )
275       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
276       CALL MPI_ALLREDUCE( sums_l(nzb,31,0), sums(nzb,31), nzt+2-nzb, &
277                           MPI_REAL, MPI_SUM, comm2d, ierr )
278       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
279       CALL MPI_ALLREDUCE( sums_l(nzb,32,0), sums(nzb,32), nzt+2-nzb, &
280                           MPI_REAL, MPI_SUM, comm2d, ierr )
281
282#else
283       sums(:,8)  = sums_l(:,8,0)
284       sums(:,30) = sums_l(:,30,0)
285       sums(:,31) = sums_l(:,31,0)
286       sums(:,32) = sums_l(:,32,0)
287#endif
288
289!
290!--    Final values are obtained by division by the total number of grid
291!--    points used for the summation.
292       hom(:,1,8,0)  = sums(:,8)  / ngp_2dh_outer(:,0)   ! e
293       hom(:,1,30,0) = sums(:,30) / ngp_2dh_outer(:,0)   ! u*2
294       hom(:,1,31,0) = sums(:,31) / ngp_2dh_outer(:,0)   ! v*2
295       hom(:,1,32,0) = sums(:,32) / ngp_2dh_outer(:,0)   ! w*2
296
297    ENDIF
298
299 END SUBROUTINE lpm_init_sgs_tke
Note: See TracBrowser for help on using the repository browser.