Home

Context Navigation

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

Last change on this file since 2016 was 2001, checked in by knoop, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 25.5 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-2016 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: diffusion_e.f90 2001 2016-08-20 18:41:22Z lvollmer $
27	!
28	! 2000 2016-08-20 18:09:15Z knoop
29	! Forced header and separation lines into 80 columns
30	!
31	! 1873 2016-04-18 14:50:06Z maronga
32	! Module renamed (removed _mod)
33	!
34	!
35	! 1850 2016-04-08 13:29:27Z maronga
36	! Module renamed
37	! Adapted for modularization of microphysics
38	!
39	! 1831 2016-04-07 13:15:51Z hoffmann
40	! turbulence renamed collision_turbulence
41	!
42	! 1682 2015-10-07 23:56:08Z knoop
43	! Code annotations made doxygen readable
44	!
45	! 1374 2014-04-25 12:55:07Z raasch
46	! missing variables added to ONLY list
47	! rif removed from acc-present-list
48	!
49	! 1340 2014-03-25 19:45:13Z kanani
50	! REAL constants defined as wp-kind
51	!
52	! 1320 2014-03-20 08:40:49Z raasch
53	! ONLY-attribute added to USE-statements,
54	! kind-parameters added to all INTEGER and REAL declaration statements,
55	! kinds are defined in new module kinds,
56	! revision history before 2012 removed,
57	! comment fields (!:) to be used for variable explanations added to
58	! all variable declaration statements
59	!
60	! 1257 2013-11-08 15:18:40Z raasch
61	! openacc loop and loop vector clauses removed
62	!
63	! 1179 2013-06-14 05:57:58Z raasch
64	! use_reference renamed use_single_reference_value
65	!
66	! 1171 2013-05-30 11:27:45Z raasch
67	! use_reference-case activated in accelerator version
68	!
69	! 1128 2013-04-12 06:19:32Z raasch
70	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
71	! j_north
72	!
73	! 1065 2012-11-22 17:42:36Z hoffmann
74	! Enabled the claculation of diss in case of turbulence = .TRUE. (parameterized
75	! effects of turbulence on autoconversion and accretion in two-moments cloud
76	! physics scheme).
77	!
78	! 1036 2012-10-22 13:43:42Z raasch
79	! code put under GPL (PALM 3.9)
80	!
81	! 1015 2012-09-27 09:23:24Z raasch
82	! accelerator version (*_acc) added,
83	! adjustment of mixing length to the Prandtl mixing length at first grid point
84	! above ground removed
85	!
86	! 1010 2012-09-20 07:59:54Z raasch
87	! cpp switch __nopointer added for pointer free version
88	!
89	! 1001 2012-09-13 14:08:46Z raasch
90	! most arrays comunicated by module instead of parameter list
91	!
92	! 825 2012-02-19 03:03:44Z raasch
93	! wang_collision_kernel renamed wang_kernel
94	!
95	! Revision 1.1 1997/09/19 07:40:24 raasch
96	! Initial revision
97	!
98	!
99	! Description:
100	! ------------
101	!> Diffusion- and dissipation terms for the TKE
102	!------------------------------------------------------------------------------!
103	MODULE diffusion_e_mod
104
105
106	PRIVATE
107	PUBLIC diffusion_e, diffusion_e_acc
108
109
110	INTERFACE diffusion_e
111	MODULE PROCEDURE diffusion_e
112	MODULE PROCEDURE diffusion_e_ij
113	END INTERFACE diffusion_e
114
115	INTERFACE diffusion_e_acc
116	MODULE PROCEDURE diffusion_e_acc
117	END INTERFACE diffusion_e_acc
118
119	CONTAINS
120
121
122	!------------------------------------------------------------------------------!
123	! Description:
124	! ------------
125	!> Call for all grid points
126	!------------------------------------------------------------------------------!
127	SUBROUTINE diffusion_e( var, var_reference )
128
129	USE arrays_3d, &
130	ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw
131
132	USE control_parameters, &
133	ONLY: atmos_ocean_sign, g, use_single_reference_value, &
134	wall_adjustment, wall_adjustment_factor
135
136	USE grid_variables, &
137	ONLY: ddx2, ddy2
138
139	USE indices, &
140	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_s_inner,&
141	nzt
142
143	USE kinds
144
145	USE microphysics_mod, &
146	ONLY: collision_turbulence
147
148	USE particle_attributes, &
149	ONLY: use_sgs_for_particles, wang_kernel
150
151	IMPLICIT NONE
152
153	INTEGER(iwp) :: i !<
154	INTEGER(iwp) :: j !<
155	INTEGER(iwp) :: k !<
156	REAL(wp) :: dvar_dz !<
157	REAL(wp) :: l_stable !<
158	REAL(wp) :: var_reference !<
159
160	#if defined( __nopointer )
161	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !<
162	#else
163	REAL(wp), DIMENSION(:,:,:), POINTER :: var !<
164	#endif
165	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dissipation !<
166	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: l !<
167	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: ll !<
168
169
170	!
171	!-- This if clause must be outside the k-loop because otherwise
172	!-- runtime errors occur with -C hopt on NEC
173	IF ( use_single_reference_value ) THEN
174
175	DO i = nxl, nxr
176	DO j = nys, nyn
177	DO k = nzb_s_inner(j,i)+1, nzt
178	!
179	!-- Calculate the mixing length (for dissipation)
180	dvar_dz = atmos_ocean_sign * &
181	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
182	IF ( dvar_dz > 0.0_wp ) THEN
183	l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
184	SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
185	ELSE
186	l_stable = l_grid(k)
187	ENDIF
188	!
189	!-- Adjustment of the mixing length
190	IF ( wall_adjustment ) THEN
191	l(k,j) = MIN( wall_adjustment_factor * &
192	( zu(k) - zw(nzb_s_inner(j,i)) ), &
193	l_grid(k), l_stable )
194	ll(k,j) = MIN( wall_adjustment_factor * &
195	( zu(k) - zw(nzb_s_inner(j,i)) ), &
196	l_grid(k) )
197	ELSE
198	l(k,j) = MIN( l_grid(k), l_stable )
199	ll(k,j) = l_grid(k)
200	ENDIF
201
202	ENDDO
203	ENDDO
204
205	!
206	!-- Calculate the tendency terms
207	DO j = nys, nyn
208	DO k = nzb_s_inner(j,i)+1, nzt
209
210	dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) * &
211	e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
212
213	tend(k,j,i) = tend(k,j,i) &
214	+ ( &
215	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
216	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
217	) * ddx2 &
218	+ ( &
219	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
220	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
221	) * ddy2 &
222	+ ( &
223	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
224	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
225	) * ddzw(k) &
226	- dissipation(k,j)
227
228	ENDDO
229	ENDDO
230
231	!
232	!-- Store dissipation if needed for calculating the sgs particle
233	!-- velocities
234	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
235	collision_turbulence ) THEN
236	DO j = nys, nyn
237	DO k = nzb_s_inner(j,i)+1, nzt
238	diss(k,j,i) = dissipation(k,j)
239	ENDDO
240	ENDDO
241	ENDIF
242
243	ENDDO
244
245	ELSE
246
247	DO i = nxl, nxr
248	DO j = nys, nyn
249	DO k = nzb_s_inner(j,i)+1, nzt
250	!
251	!-- Calculate the mixing length (for dissipation)
252	dvar_dz = atmos_ocean_sign * &
253	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
254	IF ( dvar_dz > 0.0_wp ) THEN
255	l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
256	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
257	ELSE
258	l_stable = l_grid(k)
259	ENDIF
260	!
261	!-- Adjustment of the mixing length
262	IF ( wall_adjustment ) THEN
263	l(k,j) = MIN( wall_adjustment_factor * &
264	( zu(k) - zw(nzb_s_inner(j,i)) ), &
265	l_grid(k), l_stable )
266	ll(k,j) = MIN( wall_adjustment_factor * &
267	( zu(k) - zw(nzb_s_inner(j,i)) ), &
268	l_grid(k) )
269	ELSE
270	l(k,j) = MIN( l_grid(k), l_stable )
271	ll(k,j) = l_grid(k)
272	ENDIF
273
274	ENDDO
275	ENDDO
276
277	!
278	!-- Calculate the tendency terms
279	DO j = nys, nyn
280	DO k = nzb_s_inner(j,i)+1, nzt
281
282	dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) * &
283	e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
284
285	tend(k,j,i) = tend(k,j,i) &
286	+ ( &
287	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
288	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
289	) * ddx2 &
290	+ ( &
291	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
292	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
293	) * ddy2 &
294	+ ( &
295	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
296	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
297	) * ddzw(k) &
298	- dissipation(k,j)
299
300	ENDDO
301	ENDDO
302
303	!
304	!-- Store dissipation if needed for calculating the sgs particle
305	!-- velocities
306	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
307	collision_turbulence ) THEN
308	DO j = nys, nyn
309	DO k = nzb_s_inner(j,i)+1, nzt
310	diss(k,j,i) = dissipation(k,j)
311	ENDDO
312	ENDDO
313	ENDIF
314
315	ENDDO
316
317	ENDIF
318
319	!
320	!-- Boundary condition for dissipation
321	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
322	collision_turbulence ) THEN
323	DO i = nxl, nxr
324	DO j = nys, nyn
325	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
326	ENDDO
327	ENDDO
328	ENDIF
329
330	END SUBROUTINE diffusion_e
331
332
333	!------------------------------------------------------------------------------!
334	! Description:
335	! ------------
336	!> Call for all grid points - accelerator version
337	!------------------------------------------------------------------------------!
338	SUBROUTINE diffusion_e_acc( var, var_reference )
339
340	USE arrays_3d, &
341	ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw
342
343	USE control_parameters, &
344	ONLY: atmos_ocean_sign, g, use_single_reference_value, &
345	wall_adjustment, wall_adjustment_factor
346
347	USE grid_variables, &
348	ONLY: ddx2, ddy2
349
350	USE indices, &
351	ONLY: i_left, i_right, j_north, j_south, nxlg, nxrg, nyng, nysg, &
352	nzb, nzb_s_inner, nzt
353
354	USE kinds
355
356	USE microphysics_mod, &
357	ONLY: collision_turbulence
358
359	USE particle_attributes, &
360	ONLY: use_sgs_for_particles, wang_kernel
361
362	IMPLICIT NONE
363
364	INTEGER(iwp) :: i !<
365	INTEGER(iwp) :: j !<
366	INTEGER(iwp) :: k !<
367	REAL(wp) :: dissipation !<
368	REAL(wp) :: dvar_dz !<
369	REAL(wp) :: l !<
370	REAL(wp) :: ll !<
371	REAL(wp) :: l_stable !<
372	REAL(wp) :: var_reference !<
373
374	#if defined( __nopointer )
375	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !<
376	#else
377	REAL(wp), DIMENSION(:,:,:), POINTER :: var !<
378	#endif
379
380
381	!
382	!-- This if clause must be outside the k-loop because otherwise
383	!-- runtime errors occur with -C hopt on NEC
384	IF ( use_single_reference_value ) THEN
385
386	!$acc kernels present( ddzu, ddzw, dd2zu, diss, e, km, l_grid ) &
387	!$acc present( nzb_s_inner, tend, var, zu, zw )
388	DO i = i_left, i_right
389	DO j = j_south, j_north
390	DO k = 1, nzt
391
392	IF ( k > nzb_s_inner(j,i) ) THEN
393	!
394	!-- Calculate the mixing length (for dissipation)
395	dvar_dz = atmos_ocean_sign * &
396	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
397	IF ( dvar_dz > 0.0_wp ) THEN
398	l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
399	SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
400	ELSE
401	l_stable = l_grid(k)
402	ENDIF
403	!
404	!-- Adjustment of the mixing length
405	IF ( wall_adjustment ) THEN
406	l = MIN( wall_adjustment_factor * &
407	( zu(k) - zw(nzb_s_inner(j,i)) ), &
408	l_grid(k), l_stable )
409	ll = MIN( wall_adjustment_factor * &
410	( zu(k) - zw(nzb_s_inner(j,i)) ), &
411	l_grid(k) )
412	ELSE
413	l = MIN( l_grid(k), l_stable )
414	ll = l_grid(k)
415	ENDIF
416	!
417	!-- Calculate the tendency terms
418	dissipation = ( 0.19_wp + 0.74_wp * l / ll ) * &
419	e(k,j,i) * SQRT( e(k,j,i) ) / l
420
421	tend(k,j,i) = tend(k,j,i) &
422	+ ( &
423	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
424	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
425	) * ddx2 &
426	+ ( &
427	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
428	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
429	) * ddy2 &
430	+ ( &
431	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
432	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
433	) * ddzw(k) &
434	- dissipation
435
436	!
437	!-- Store dissipation if needed for calculating the sgs particle
438	!-- velocities
439	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
440	collision_turbulence ) THEN
441	diss(k,j,i) = dissipation
442	ENDIF
443
444	ENDIF
445
446	ENDDO
447	ENDDO
448	ENDDO
449	!$acc end kernels
450
451	ELSE
452
453	!$acc kernels present( ddzu, ddzw, dd2zu, diss, e, km, l_grid ) &
454	!$acc present( nzb_s_inner, tend, var, zu, zw )
455	DO i = i_left, i_right
456	DO j = j_south, j_north
457	DO k = 1, nzt
458
459	IF ( k > nzb_s_inner(j,i) ) THEN
460	!
461	!-- Calculate the mixing length (for dissipation)
462	dvar_dz = atmos_ocean_sign * &
463	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
464	IF ( dvar_dz > 0.0_wp ) THEN
465	l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
466	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
467	ELSE
468	l_stable = l_grid(k)
469	ENDIF
470	!
471	!-- Adjustment of the mixing length
472	IF ( wall_adjustment ) THEN
473	l = MIN( wall_adjustment_factor * &
474	( zu(k) - zw(nzb_s_inner(j,i)) ), &
475	l_grid(k), l_stable )
476	ll = MIN( wall_adjustment_factor * &
477	( zu(k) - zw(nzb_s_inner(j,i)) ), &
478	l_grid(k) )
479	ELSE
480	l = MIN( l_grid(k), l_stable )
481	ll = l_grid(k)
482	ENDIF
483	!
484	!-- Calculate the tendency terms
485	dissipation = ( 0.19_wp + 0.74_wp * l / ll ) * &
486	e(k,j,i) * SQRT( e(k,j,i) ) / l
487
488	tend(k,j,i) = tend(k,j,i) &
489	+ ( &
490	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
491	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
492	) * ddx2 &
493	+ ( &
494	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
495	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
496	) * ddy2 &
497	+ ( &
498	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
499	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
500	) * ddzw(k) &
501	- dissipation
502
503	!
504	!-- Store dissipation if needed for calculating the sgs
505	!-- particle velocities
506	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
507	collision_turbulence ) THEN
508	diss(k,j,i) = dissipation
509	ENDIF
510
511	ENDIF
512
513	ENDDO
514	ENDDO
515	ENDDO
516	!$acc end kernels
517
518	ENDIF
519
520	!
521	!-- Boundary condition for dissipation
522	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
523	collision_turbulence ) THEN
524	!$acc kernels present( diss, nzb_s_inner )
525	DO i = i_left, i_right
526	DO j = j_south, j_north
527	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
528	ENDDO
529	ENDDO
530	!$acc end kernels
531	ENDIF
532
533	END SUBROUTINE diffusion_e_acc
534
535
536	!------------------------------------------------------------------------------!
537	! Description:
538	! ------------
539	!> Call for grid point i,j
540	!------------------------------------------------------------------------------!
541	SUBROUTINE diffusion_e_ij( i, j, var, var_reference )
542
543	USE arrays_3d, &
544	ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw
545
546	USE control_parameters, &
547	ONLY: atmos_ocean_sign, g, use_single_reference_value, &
548	wall_adjustment, wall_adjustment_factor
549
550	USE grid_variables, &
551	ONLY: ddx2, ddy2
552
553	USE indices, &
554	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzb_s_inner, nzt
555
556	USE kinds
557
558	USE microphysics_mod, &
559	ONLY: collision_turbulence
560
561	USE particle_attributes, &
562	ONLY: use_sgs_for_particles, wang_kernel
563
564	IMPLICIT NONE
565
566	INTEGER(iwp) :: i !<
567	INTEGER(iwp) :: j !<
568	INTEGER(iwp) :: k !<
569	REAL(wp) :: dvar_dz !<
570	REAL(wp) :: l_stable !<
571	REAL(wp) :: var_reference !<
572
573	#if defined( __nopointer )
574	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !<
575	#else
576	REAL(wp), DIMENSION(:,:,:), POINTER :: var !<
577	#endif
578	REAL(wp), DIMENSION(nzb+1:nzt) :: dissipation !<
579	REAL(wp), DIMENSION(nzb+1:nzt) :: l !<
580	REAL(wp), DIMENSION(nzb+1:nzt) :: ll !<
581
582
583	!
584	!-- Calculate the mixing length (for dissipation)
585	DO k = nzb_s_inner(j,i)+1, nzt
586	dvar_dz = atmos_ocean_sign * &
587	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
588	IF ( dvar_dz > 0.0_wp ) THEN
589	IF ( use_single_reference_value ) THEN
590	l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
591	SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
592	ELSE
593	l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
594	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
595	ENDIF
596	ELSE
597	l_stable = l_grid(k)
598	ENDIF
599	!
600	!-- Adjustment of the mixing length
601	IF ( wall_adjustment ) THEN
602	l(k) = MIN( wall_adjustment_factor * &
603	( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k), &
604	l_stable )
605	ll(k) = MIN( wall_adjustment_factor * &
606	( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k) )
607	ELSE
608	l(k) = MIN( l_grid(k), l_stable )
609	ll(k) = l_grid(k)
610	ENDIF
611	!
612	!-- Calculate the tendency term
613	dissipation(k) = ( 0.19_wp + 0.74_wp * l(k) / ll(k) ) * e(k,j,i) * &
614	SQRT( e(k,j,i) ) / l(k)
615
616	tend(k,j,i) = tend(k,j,i) &
617	+ ( &
618	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
619	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
620	) * ddx2 &
621	+ ( &
622	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
623	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
624	) * ddy2 &
625	+ ( &
626	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
627	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
628	) * ddzw(k) &
629	- dissipation(k)
630
631	ENDDO
632
633	!
634	!-- Store dissipation if needed for calculating the sgs particle velocities
635	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
636	collision_turbulence ) THEN
637	DO k = nzb_s_inner(j,i)+1, nzt
638	diss(k,j,i) = dissipation(k)
639	ENDDO
640	!
641	!-- Boundary condition for dissipation
642	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
643	ENDIF
644
645	END SUBROUTINE diffusion_e_ij
646
647	END MODULE diffusion_e_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |