Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |