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

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

Rephrase comment

  • Property svn:keywords set to Id
File size: 12.9 KB
Line 
1!> @file lpm_init_sgs_tke.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
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-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: lpm_init_sgs_tke.f90 3985 2019-05-16 17:00:37Z suehring $
27! Improve comment
28!
29! 3981 2019-05-15 14:43:01Z suehring
30! Bugfix in particle nesting, set boundary conditions at non-cyclic boundaries
31!
32! 3655 2019-01-07 16:51:22Z knoop
33! unused variables removed
34!
35! 2718 2018-01-02 08:49:38Z maronga
36! Corrected "Former revisions" section
37!
38! 2696 2017-12-14 17:12:51Z kanani
39! Change in file header (GPL part)
40!
41! 2233 2017-05-30 18:08:54Z suehring
42!
43! 2232 2017-05-30 17:47:52Z suehring
44! Adjustments according to new topography realization
45!
46! 2000 2016-08-20 18:09:15Z knoop
47! Forced header and separation lines into 80 columns
48!
49! 1929 2016-06-09 16:25:25Z suehring
50! sgs_wfu_par, sgs_wfv_par and sgs_wfw_par are replaced by sgs_wf_par
51!
52! 1822 2016-04-07 07:49:42Z hoffmann
53! Unused variables removed.
54!
55! 1682 2015-10-07 23:56:08Z knoop
56! Code annotations made doxygen readable
57!
58! 1359 2014-04-11 17:15:14Z hoffmann
59! New particle structure integrated.
60! Kind definition added to all floating point numbers.
61!
62! 1320 2014-03-20 08:40:49Z raasch
63! ONLY-attribute added to USE-statements,
64! kind-parameters added to all INTEGER and REAL declaration statements,
65! kinds are defined in new module kinds,
66! comment fields (!:) to be used for variable explanations added to
67! all variable declaration statements
68!
69! 1036 2012-10-22 13:43:42Z raasch
70! code put under GPL (PALM 3.9)
71!
72! 849 2012-03-15 10:35:09Z raasch
73! initial revision (former part of advec_particles)
74!
75!
76! Description:
77! ------------
78!> Calculates quantities required for considering the SGS velocity fluctuations
79!> in the particle transport by a stochastic approach. The respective
80!> quantities are: SGS-TKE gradients and horizontally averaged profiles of the
81!> SGS TKE and the resolved-scale velocity variances.
82!------------------------------------------------------------------------------!
83 SUBROUTINE lpm_init_sgs_tke
84 
85
86    USE arrays_3d,                                                             &
87        ONLY:  de_dx, de_dy, de_dz, diss, e, u, v, w, zu
88
89    USE control_parameters,                                                    &
90        ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s
91       
92    USE grid_variables,                                                        &
93        ONLY:  ddx, ddy
94
95    USE indices,                                                               &
96        ONLY:  nbgp, ngp_2dh_outer, nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
97
98    USE kinds
99
100    USE particle_attributes,                                                   &
101        ONLY:  sgs_wf_part
102
103    USE pegrid
104
105    USE statistics,                                                            &
106        ONLY:  flow_statistics_called, hom, sums, sums_l
107
108    USE surface_mod,                                                           &
109        ONLY:  bc_h
110
111    IMPLICIT NONE
112
113    INTEGER(iwp) ::  i      !< index variable along x
114    INTEGER(iwp) ::  j      !< index variable along y
115    INTEGER(iwp) ::  k      !< index variable along z
116    INTEGER(iwp) ::  m      !< running index for the surface elements
117
118    REAL(wp) ::  flag1      !< flag to mask topography
119
120!
121!-- TKE gradient along x and y
122    DO  i = nxl, nxr
123       DO  j = nys, nyn
124          DO  k = nzb, nzt+1
125
126             IF ( .NOT. BTEST( wall_flags_0(k,j,i-1), 0 )  .AND.               &
127                        BTEST( wall_flags_0(k,j,i), 0   )  .AND.               &
128                        BTEST( wall_flags_0(k,j,i+1), 0 ) )                    &
129             THEN
130                de_dx(k,j,i) = 2.0_wp * sgs_wf_part *                          &
131                               ( e(k,j,i+1) - e(k,j,i) ) * ddx
132             ELSEIF ( BTEST( wall_flags_0(k,j,i-1), 0 )  .AND.                 &
133                      BTEST( wall_flags_0(k,j,i), 0   )  .AND.                 &
134                .NOT. BTEST( wall_flags_0(k,j,i+1), 0 ) )                      &
135             THEN
136                de_dx(k,j,i) = 2.0_wp * sgs_wf_part *                          &
137                               ( e(k,j,i) - e(k,j,i-1) ) * ddx
138             ELSEIF ( .NOT. BTEST( wall_flags_0(k,j,i), 22   )  .AND.          &
139                      .NOT. BTEST( wall_flags_0(k,j,i+1), 22 ) )               &   
140             THEN
141                de_dx(k,j,i) = 0.0_wp
142             ELSEIF ( .NOT. BTEST( wall_flags_0(k,j,i-1), 22 )  .AND.          &
143                      .NOT. BTEST( wall_flags_0(k,j,i), 22   ) )               &
144             THEN
145                de_dx(k,j,i) = 0.0_wp
146             ELSE
147                de_dx(k,j,i) = sgs_wf_part * ( e(k,j,i+1) - e(k,j,i-1) ) * ddx
148             ENDIF
149
150             IF ( .NOT. BTEST( wall_flags_0(k,j-1,i), 0 )  .AND.               &
151                        BTEST( wall_flags_0(k,j,i), 0   )  .AND.               &
152                        BTEST( wall_flags_0(k,j+1,i), 0 ) )                    &
153             THEN
154                de_dy(k,j,i) = 2.0_wp * sgs_wf_part *                          &
155                               ( e(k,j+1,i) - e(k,j,i) ) * ddy
156             ELSEIF ( BTEST( wall_flags_0(k,j-1,i), 0 )  .AND.                 &
157                      BTEST( wall_flags_0(k,j,i), 0   )  .AND.                 &
158                .NOT. BTEST( wall_flags_0(k,j+1,i), 0 ) )                      &
159             THEN
160                de_dy(k,j,i) = 2.0_wp * sgs_wf_part *                          &
161                               ( e(k,j,i) - e(k,j-1,i) ) * ddy
162             ELSEIF ( .NOT. BTEST( wall_flags_0(k,j,i), 22   )  .AND.          &
163                      .NOT. BTEST( wall_flags_0(k,j+1,i), 22 ) )               &   
164             THEN
165                de_dy(k,j,i) = 0.0_wp
166             ELSEIF ( .NOT. BTEST( wall_flags_0(k,j-1,i), 22 )  .AND.          &
167                      .NOT. BTEST( wall_flags_0(k,j,i), 22   ) )               &
168             THEN
169                de_dy(k,j,i) = 0.0_wp
170             ELSE
171                de_dy(k,j,i) = sgs_wf_part * ( e(k,j+1,i) - e(k,j-1,i) ) * ddy
172             ENDIF
173
174          ENDDO
175       ENDDO
176    ENDDO
177
178!
179!-- TKE gradient along z at topograhy and  including bottom and top boundary conditions
180    DO  i = nxl, nxr
181       DO  j = nys, nyn
182          DO  k = nzb+1, nzt-1
183!
184!--          Flag to mask topography
185             flag1 = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0  ) )
186
187             de_dz(k,j,i) = 2.0_wp * sgs_wf_part *                             &
188                           ( e(k+1,j,i) - e(k-1,j,i) ) / ( zu(k+1) - zu(k-1) ) &
189                                                 * flag1 
190          ENDDO
191!
192!--       upward-facing surfaces
193          DO  m = bc_h(0)%start_index(j,i), bc_h(0)%end_index(j,i)
194             k            = bc_h(0)%k(m)
195             de_dz(k,j,i) = 2.0_wp * sgs_wf_part *                             &
196                           ( e(k+1,j,i) - e(k,j,i)   ) / ( zu(k+1) - zu(k) )
197          ENDDO
198!
199!--       downward-facing surfaces
200          DO  m = bc_h(1)%start_index(j,i), bc_h(1)%end_index(j,i)
201             k            = bc_h(1)%k(m)
202             de_dz(k,j,i) = 2.0_wp * sgs_wf_part *                             &
203                           ( e(k,j,i) - e(k-1,j,i)   ) / ( zu(k) - zu(k-1) )
204          ENDDO
205
206          de_dz(nzb,j,i)   = 0.0_wp
207          de_dz(nzt,j,i)   = 0.0_wp
208          de_dz(nzt+1,j,i) = 0.0_wp
209       ENDDO
210    ENDDO
211!
212!-- Ghost point exchange
213    CALL exchange_horiz( de_dx, nbgp )
214    CALL exchange_horiz( de_dy, nbgp )
215    CALL exchange_horiz( de_dz, nbgp )
216    CALL exchange_horiz( diss, nbgp  )
217!
218!-- Set boundary conditions at non-periodic boundaries. Note, at non-period
219!-- boundaries zero-gradient boundary conditions are set for the subgrid TKE.
220!-- Thus, TKE gradients normal to the respective lateral boundaries are zero,
221!-- while tangetial TKE gradients then must be the same as within the prognostic
222!-- domain. 
223    IF ( bc_dirichlet_l )  THEN
224       de_dx(:,:,-1) = 0.0_wp
225       de_dy(:,:,-1) = de_dy(:,:,0) 
226       de_dz(:,:,-1) = de_dz(:,:,0)
227    ENDIF
228    IF ( bc_dirichlet_r )  THEN
229       de_dx(:,:,nxr+1) = 0.0_wp
230       de_dy(:,:,nxr+1) = de_dy(:,:,nxr) 
231       de_dz(:,:,nxr+1) = de_dz(:,:,nxr)
232    ENDIF
233    IF ( bc_dirichlet_n )  THEN
234       de_dx(:,nyn+1,:) = de_dx(:,nyn,:)
235       de_dy(:,nyn+1,:) = 0.0_wp 
236       de_dz(:,nyn+1,:) = de_dz(:,nyn,:)
237    ENDIF
238    IF ( bc_dirichlet_s )  THEN
239       de_dx(:,nys-1,:) = de_dx(:,nys,:)
240       de_dy(:,nys-1,:) = 0.0_wp 
241       de_dz(:,nys-1,:) = de_dz(:,nys,:)
242    ENDIF 
243!
244!-- Calculate the horizontally averaged profiles of SGS TKE and resolved
245!-- velocity variances (they may have been already calculated in routine
246!-- flow_statistics).
247    IF ( .NOT. flow_statistics_called )  THEN
248
249!
250!--    First calculate horizontally averaged profiles of the horizontal
251!--    velocities.
252       sums_l(:,1,0) = 0.0_wp
253       sums_l(:,2,0) = 0.0_wp
254
255       DO  i = nxl, nxr
256          DO  j =  nys, nyn
257             DO  k = nzb, nzt+1
258!
259!--             Flag indicating vicinity of wall
260                flag1 = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 24 ) )
261
262                sums_l(k,1,0)  = sums_l(k,1,0)  + u(k,j,i) * flag1
263                sums_l(k,2,0)  = sums_l(k,2,0)  + v(k,j,i) * flag1
264             ENDDO
265          ENDDO
266       ENDDO
267
268#if defined( __parallel )
269!
270!--    Compute total sum from local sums
271       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
272       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, &
273                           MPI_REAL, MPI_SUM, comm2d, ierr )
274       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
275       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, &
276                              MPI_REAL, MPI_SUM, comm2d, ierr )
277#else
278       sums(:,1) = sums_l(:,1,0)
279       sums(:,2) = sums_l(:,2,0)
280#endif
281
282!
283!--    Final values are obtained by division by the total number of grid
284!--    points used for the summation.
285       hom(:,1,1,0) = sums(:,1) / ngp_2dh_outer(:,0)   ! u
286       hom(:,1,2,0) = sums(:,2) / ngp_2dh_outer(:,0)   ! v
287
288!
289!--    Now calculate the profiles of SGS TKE and the resolved-scale
290!--    velocity variances
291       sums_l(:,8,0)  = 0.0_wp
292       sums_l(:,30,0) = 0.0_wp
293       sums_l(:,31,0) = 0.0_wp
294       sums_l(:,32,0) = 0.0_wp
295       DO  i = nxl, nxr
296          DO  j = nys, nyn
297             DO  k = nzb, nzt+1
298!
299!--             Flag indicating vicinity of wall
300                flag1 = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 24 ) )
301
302                sums_l(k,8,0)  = sums_l(k,8,0)  + e(k,j,i)                       * flag1
303                sums_l(k,30,0) = sums_l(k,30,0) + ( u(k,j,i) - hom(k,1,1,0) )**2 * flag1
304                sums_l(k,31,0) = sums_l(k,31,0) + ( v(k,j,i) - hom(k,1,2,0) )**2 * flag1
305                sums_l(k,32,0) = sums_l(k,32,0) + w(k,j,i)**2                    * flag1
306             ENDDO
307          ENDDO
308       ENDDO
309
310#if defined( __parallel )
311!
312!--    Compute total sum from local sums
313       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
314       CALL MPI_ALLREDUCE( sums_l(nzb,8,0), sums(nzb,8), nzt+2-nzb, &
315                           MPI_REAL, MPI_SUM, comm2d, ierr )
316       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
317       CALL MPI_ALLREDUCE( sums_l(nzb,30,0), sums(nzb,30), nzt+2-nzb, &
318                           MPI_REAL, MPI_SUM, comm2d, ierr )
319       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
320       CALL MPI_ALLREDUCE( sums_l(nzb,31,0), sums(nzb,31), nzt+2-nzb, &
321                           MPI_REAL, MPI_SUM, comm2d, ierr )
322       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
323       CALL MPI_ALLREDUCE( sums_l(nzb,32,0), sums(nzb,32), nzt+2-nzb, &
324                           MPI_REAL, MPI_SUM, comm2d, ierr )
325
326#else
327       sums(:,8)  = sums_l(:,8,0)
328       sums(:,30) = sums_l(:,30,0)
329       sums(:,31) = sums_l(:,31,0)
330       sums(:,32) = sums_l(:,32,0)
331#endif
332
333!
334!--    Final values are obtained by division by the total number of grid
335!--    points used for the summation.
336       hom(:,1,8,0)  = sums(:,8)  / ngp_2dh_outer(:,0)   ! e
337       hom(:,1,30,0) = sums(:,30) / ngp_2dh_outer(:,0)   ! u*2
338       hom(:,1,31,0) = sums(:,31) / ngp_2dh_outer(:,0)   ! v*2
339       hom(:,1,32,0) = sums(:,32) / ngp_2dh_outer(:,0)   ! w*2
340
341    ENDIF
342
343 END SUBROUTINE lpm_init_sgs_tke
Note: See TracBrowser for help on using the repository browser.