Home

Context Navigation

source: palm/trunk/SOURCE/diffusion_e.f90 @ 2605

Last change on this file since 2605 was 2233, checked in by suehring, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 19.7 KB

Line
1	!> @file diffusion_e.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	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: diffusion_e.f90 2233 2017-05-30 18:08:54Z raasch $
27	!
28	! 2232 2017-05-30 17:47:52Z suehring
29	! Adjustments to new topography and surface concept
30	!
31	! 2118 2017-01-17 16:38:49Z raasch
32	! OpenACC version of subroutine removed
33	!
34	! 2037 2016-10-26 11:15:40Z knoop
35	! Anelastic approximation implemented
36	!
37	! 2000 2016-08-20 18:09:15Z knoop
38	! Forced header and separation lines into 80 columns
39	!
40	! 1873 2016-04-18 14:50:06Z maronga
41	! Module renamed (removed _mod)
42	!
43	! 1850 2016-04-08 13:29:27Z maronga
44	! Module renamed
45	! Adapted for modularization of microphysics
46	!
47	! 1831 2016-04-07 13:15:51Z hoffmann
48	! turbulence renamed collision_turbulence
49	!
50	! 1682 2015-10-07 23:56:08Z knoop
51	! Code annotations made doxygen readable
52	!
53	! 1374 2014-04-25 12:55:07Z raasch
54	! missing variables added to ONLY list
55	! rif removed from acc-present-list
56	!
57	! 1340 2014-03-25 19:45:13Z kanani
58	! REAL constants defined as wp-kind
59	!
60	! 1320 2014-03-20 08:40:49Z raasch
61	! ONLY-attribute added to USE-statements,
62	! kind-parameters added to all INTEGER and REAL declaration statements,
63	! kinds are defined in new module kinds,
64	! revision history before 2012 removed,
65	! comment fields (!:) to be used for variable explanations added to
66	! all variable declaration statements
67	!
68	! 1257 2013-11-08 15:18:40Z raasch
69	! openacc loop and loop vector clauses removed
70	!
71	! 1179 2013-06-14 05:57:58Z raasch
72	! use_reference renamed use_single_reference_value
73	!
74	! 1171 2013-05-30 11:27:45Z raasch
75	! use_reference-case activated in accelerator version
76	!
77	! 1128 2013-04-12 06:19:32Z raasch
78	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
79	! j_north
80	!
81	! 1065 2012-11-22 17:42:36Z hoffmann
82	! Enabled the claculation of diss in case of turbulence = .TRUE. (parameterized
83	! effects of turbulence on autoconversion and accretion in two-moments cloud
84	! physics scheme).
85	!
86	! 1036 2012-10-22 13:43:42Z raasch
87	! code put under GPL (PALM 3.9)
88	!
89	! 1015 2012-09-27 09:23:24Z raasch
90	! accelerator version (*_acc) added,
91	! adjustment of mixing length to the Prandtl mixing length at first grid point
92	! above ground removed
93	!
94	! 1010 2012-09-20 07:59:54Z raasch
95	! cpp switch __nopointer added for pointer free version
96	!
97	! 1001 2012-09-13 14:08:46Z raasch
98	! most arrays comunicated by module instead of parameter list
99	!
100	! 825 2012-02-19 03:03:44Z raasch
101	! wang_collision_kernel renamed wang_kernel
102	!
103	! Revision 1.1 1997/09/19 07:40:24 raasch
104	! Initial revision
105	!
106	!
107	! Description:
108	! ------------
109	!> Diffusion- and dissipation terms for the TKE
110	!------------------------------------------------------------------------------!
111	MODULE diffusion_e_mod
112
113
114	PRIVATE
115	PUBLIC diffusion_e
116
117
118	INTERFACE diffusion_e
119	MODULE PROCEDURE diffusion_e
120	MODULE PROCEDURE diffusion_e_ij
121	END INTERFACE diffusion_e
122
123	CONTAINS
124
125
126	!------------------------------------------------------------------------------!
127	! Description:
128	! ------------
129	!> Call for all grid points
130	!------------------------------------------------------------------------------!
131	SUBROUTINE diffusion_e( var, var_reference )
132
133	USE arrays_3d, &
134	ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, l_wall, tend, zu, zw,&
135	drho_air, rho_air_zw
136
137	USE control_parameters, &
138	ONLY: atmos_ocean_sign, g, use_single_reference_value, &
139	wall_adjustment, wall_adjustment_factor
140
141	USE grid_variables, &
142	ONLY: ddx2, ddy2
143
144	USE indices, &
145	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_max, &
146	nzt, wall_flags_0
147
148	USE kinds
149
150	USE microphysics_mod, &
151	ONLY: collision_turbulence
152
153	USE particle_attributes, &
154	ONLY: use_sgs_for_particles, wang_kernel
155
156	USE surface_mod, &
157	ONLY : bc_h
158
159	IMPLICIT NONE
160
161	INTEGER(iwp) :: i !< running index x direction
162	INTEGER(iwp) :: j !< running index y direction
163	INTEGER(iwp) :: k !< running index z direction
164	INTEGER(iwp) :: m !< running index surface elements
165
166	REAL(wp) :: dvar_dz !<
167	REAL(wp) :: flag !< flag to mask topography
168	REAL(wp) :: l_stable !<
169	REAL(wp) :: var_reference !<
170
171	#if defined( __nopointer )
172	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !<
173	#else
174	REAL(wp), DIMENSION(:,:,:), POINTER :: var !<
175	#endif
176	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dissipation !<
177	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: l !<
178	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: ll !<
179
180
181	!
182	!-- This if clause must be outside the k-loop because otherwise
183	!-- runtime errors occur with -C hopt on NEC
184	IF ( use_single_reference_value ) THEN
185
186	DO i = nxl, nxr
187	DO j = nys, nyn
188	DO k = nzb+1, nzt
189	!
190	!-- Calculate the mixing length (for dissipation)
191	dvar_dz = atmos_ocean_sign * &
192	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
193	IF ( dvar_dz > 0.0_wp ) THEN
194	l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
195	SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
196	ELSE
197	l_stable = l_grid(k)
198	ENDIF
199	!
200	!-- Adjustment of the mixing length
201	IF ( wall_adjustment ) THEN
202	l(k,j) = MIN( wall_adjustment_factor * l_wall(k,j,i), &
203	l_grid(k), l_stable )
204	ll(k,j) = MIN( wall_adjustment_factor * l_wall(k,j,i), &
205	l_grid(k) )
206	ELSE
207	l(k,j) = MIN( l_grid(k), l_stable )
208	ll(k,j) = l_grid(k)
209	ENDIF
210
211	ENDDO
212	ENDDO
213
214	!
215	!-- Calculate the tendency terms
216	DO j = nys, nyn
217	DO k = nzb+1, nzt
218	!
219	!-- Predetermine flag to mask topography
220	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
221
222	dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) * &
223	e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
224
225	tend(k,j,i) = tend(k,j,i) &
226	+ ( &
227	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
228	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
229	) * ddx2 * flag &
230	+ ( &
231	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
232	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
233	) * ddy2 * flag &
234	+ ( &
235	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
236	* rho_air_zw(k) &
237	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
238	* rho_air_zw(k-1) &
239	) * ddzw(k) * drho_air(k) * flag &
240	- dissipation(k,j) * flag
241
242	ENDDO
243	ENDDO
244
245	!
246	!-- Store dissipation if needed for calculating the sgs particle
247	!-- velocities
248	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
249	collision_turbulence ) THEN
250	DO j = nys, nyn
251	DO k = nzb+1, nzt
252	diss(k,j,i) = dissipation(k,j) * &
253	MERGE( 1.0_wp, 0.0_wp, &
254	BTEST( wall_flags_0(k,j,i), 0 ) )
255	ENDDO
256	ENDDO
257	ENDIF
258
259	ENDDO
260
261	ELSE
262
263	DO i = nxl, nxr
264	DO j = nys, nyn
265	DO k = nzb+1, nzt
266	!
267	!-- Calculate the mixing length (for dissipation)
268	dvar_dz = atmos_ocean_sign * &
269	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
270	IF ( dvar_dz > 0.0_wp ) THEN
271	l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
272	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
273	ELSE
274	l_stable = l_grid(k)
275	ENDIF
276	!
277	!-- Adjustment of the mixing length
278	IF ( wall_adjustment ) THEN
279	l(k,j) = MIN( wall_adjustment_factor * l_wall(k,j,i), &
280	l_grid(k), l_stable )
281	ll(k,j) = MIN( wall_adjustment_factor * l_wall(k,j,i), &
282	l_grid(k) )
283	ELSE
284	l(k,j) = MIN( l_grid(k), l_stable )
285	ll(k,j) = l_grid(k)
286	ENDIF
287
288	ENDDO
289	ENDDO
290
291	!
292	!-- Calculate the tendency terms
293	DO j = nys, nyn
294	DO k = nzb+1, nzt
295	!
296	!-- Predetermine flag to mask topography
297	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
298
299	dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) *&
300	e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
301
302	tend(k,j,i) = tend(k,j,i) &
303	+ ( &
304	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
305	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
306	) * ddx2 * flag &
307	+ ( &
308	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
309	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
310	) * ddy2 * flag &
311	+ ( &
312	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
313	* rho_air_zw(k) &
314	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
315	* rho_air_zw(k-1) &
316	) * ddzw(k) * drho_air(k) * flag &
317	- dissipation(k,j) * flag
318
319	ENDDO
320	ENDDO
321
322	!
323	!-- Store dissipation if needed for calculating the sgs particle
324	!-- velocities
325	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
326	collision_turbulence ) THEN
327	DO j = nys, nyn
328	DO k = nzb+1, nzt
329	diss(k,j,i) = dissipation(k,j) * &
330	MERGE( 1.0_wp, 0.0_wp, &
331	BTEST( wall_flags_0(k,j,i), 0 ) )
332	ENDDO
333	ENDDO
334	ENDIF
335
336	ENDDO
337
338	ENDIF
339
340	!
341	!-- Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:)
342	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
343	collision_turbulence ) THEN
344	!
345	!-- Upward facing surfaces
346	DO m = 1, bc_h(0)%ns
347	i = bc_h(0)%i(m)
348	j = bc_h(0)%j(m)
349	k = bc_h(0)%k(m)
350	diss(k-1,j,i) = diss(k,j,i)
351	ENDDO
352	!
353	!-- Downward facing surfaces
354	DO m = 1, bc_h(1)%ns
355	i = bc_h(1)%i(m)
356	j = bc_h(1)%j(m)
357	k = bc_h(1)%k(m)
358	diss(k+1,j,i) = diss(k,j,i)
359	ENDDO
360
361	ENDIF
362
363	END SUBROUTINE diffusion_e
364
365
366	!------------------------------------------------------------------------------!
367	! Description:
368	! ------------
369	!> Call for grid point i,j
370	!------------------------------------------------------------------------------!
371	SUBROUTINE diffusion_e_ij( i, j, var, var_reference )
372
373	USE arrays_3d, &
374	ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, l_wall, tend, zu, zw,&
375	drho_air, rho_air_zw
376
377	USE control_parameters, &
378	ONLY: atmos_ocean_sign, g, use_single_reference_value, &
379	wall_adjustment, wall_adjustment_factor
380
381	USE grid_variables, &
382	ONLY: ddx2, ddy2
383
384	USE indices, &
385	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzb_max, nzt, wall_flags_0
386
387	USE kinds
388
389	USE microphysics_mod, &
390	ONLY: collision_turbulence
391
392	USE particle_attributes, &
393	ONLY: use_sgs_for_particles, wang_kernel
394
395	USE surface_mod, &
396	ONLY : bc_h
397
398	IMPLICIT NONE
399
400	INTEGER(iwp) :: i !< running index x direction
401	INTEGER(iwp) :: j !< running index y direction
402	INTEGER(iwp) :: k !< running index z direction
403	INTEGER(iwp) :: m !< running index surface elements
404	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
405	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
406
407	REAL(wp) :: dvar_dz !<
408	REAL(wp) :: flag !< flag to mask topography
409	REAL(wp) :: l_stable !<
410	REAL(wp) :: var_reference !<
411
412	#if defined( __nopointer )
413	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !<
414	#else
415	REAL(wp), DIMENSION(:,:,:), POINTER :: var !<
416	#endif
417	REAL(wp), DIMENSION(nzb+1:nzt) :: dissipation !<
418	REAL(wp), DIMENSION(nzb+1:nzt) :: l !<
419	REAL(wp), DIMENSION(nzb+1:nzt) :: ll !<
420
421	!
422	!-- Calculate the mixing length (for dissipation)
423	DO k = nzb+1, nzt
424	!
425	!-- Predetermine flag to mask topography
426	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
427
428	dvar_dz = atmos_ocean_sign * &
429	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
430	IF ( dvar_dz > 0.0_wp ) THEN
431	IF ( use_single_reference_value ) THEN
432	l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
433	SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
434	ELSE
435	l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
436	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
437	ENDIF
438	ELSE
439	l_stable = l_grid(k)
440	ENDIF
441	!
442	!-- Adjustment of the mixing length
443	IF ( wall_adjustment ) THEN
444	l(k) = MIN( wall_adjustment_factor * l_wall(k,j,i), &
445	l_grid(k), l_stable )
446	ll(k) = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k) )
447	ELSE
448	l(k) = MIN( l_grid(k), l_stable )
449	ll(k) = l_grid(k)
450	ENDIF
451	!
452	!-- Calculate the tendency term
453	dissipation(k) = ( 0.19_wp + 0.74_wp * l(k) / ll(k) ) * e(k,j,i) * &
454	SQRT( e(k,j,i) ) / l(k)
455
456	tend(k,j,i) = tend(k,j,i) &
457	+ ( &
458	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
459	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
460	) * ddx2 * flag &
461	+ ( &
462	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
463	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
464	) * ddy2 * flag &
465	+ ( &
466	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
467	* rho_air_zw(k) &
468	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
469	* rho_air_zw(k-1) &
470	) * ddzw(k) * drho_air(k) * flag &
471	- dissipation(k) * flag
472
473	ENDDO
474
475	!
476	!-- Store dissipation if needed for calculating the sgs particle velocities
477	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
478	collision_turbulence ) THEN
479	DO k = nzb+1, nzt
480	diss(k,j,i) = dissipation(k) * &
481	MERGE( 1.0_wp, 0.0_wp, &
482	BTEST( wall_flags_0(k,j,i), 0 ) )
483	ENDDO
484	!
485	!-- Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:)
486	!-- For each surface type determine start and end index (in case of elevated
487	!-- toopography several up/downward facing surfaces may exist.
488	surf_s = bc_h(0)%start_index(j,i)
489	surf_e = bc_h(0)%end_index(j,i)
490	DO m = surf_s, surf_e
491	k = bc_h(0)%k(m)
492	diss(k-1,j,i) = diss(k,j,i)
493	ENDDO
494	!
495	!-- Downward facing surfaces
496	surf_s = bc_h(1)%start_index(j,i)
497	surf_e = bc_h(1)%end_index(j,i)
498	DO m = surf_s, surf_e
499	k = bc_h(1)%k(m)
500	diss(k+1,j,i) = diss(k,j,i)
501	ENDDO
502	ENDIF
503
504	END SUBROUTINE diffusion_e_ij
505
506	END MODULE diffusion_e_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |