Home

Context Navigation

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

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |