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

Last change on this file since 3117 was 2718, checked in by maronga, 7 years ago

deleting of deprecated files; headers updated where needed

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