Home

Context Navigation

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

Last change on this file since 2037 was 2037, checked in by knoop, 8 years ago
Anelastic approximation implemented
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	! Anelastic approximation implemented
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: diffusion_e.f90 2037 2016-10-26 11:15:40Z knoop $
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	drho_air, rho_air_zw
132
133	USE control_parameters, &
134	ONLY: atmos_ocean_sign, g, use_single_reference_value, &
135	wall_adjustment, wall_adjustment_factor
136
137	USE grid_variables, &
138	ONLY: ddx2, ddy2
139
140	USE indices, &
141	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_s_inner,&
142	nzt
143
144	USE kinds
145
146	USE microphysics_mod, &
147	ONLY: collision_turbulence
148
149	USE particle_attributes, &
150	ONLY: use_sgs_for_particles, wang_kernel
151
152	IMPLICIT NONE
153
154	INTEGER(iwp) :: i !<
155	INTEGER(iwp) :: j !<
156	INTEGER(iwp) :: k !<
157	REAL(wp) :: dvar_dz !<
158	REAL(wp) :: l_stable !<
159	REAL(wp) :: var_reference !<
160
161	#if defined( __nopointer )
162	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !<
163	#else
164	REAL(wp), DIMENSION(:,:,:), POINTER :: var !<
165	#endif
166	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dissipation !<
167	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: l !<
168	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: ll !<
169
170
171	!
172	!-- This if clause must be outside the k-loop because otherwise
173	!-- runtime errors occur with -C hopt on NEC
174	IF ( use_single_reference_value ) THEN
175
176	DO i = nxl, nxr
177	DO j = nys, nyn
178	DO k = nzb_s_inner(j,i)+1, nzt
179	!
180	!-- Calculate the mixing length (for dissipation)
181	dvar_dz = atmos_ocean_sign * &
182	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
183	IF ( dvar_dz > 0.0_wp ) THEN
184	l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
185	SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
186	ELSE
187	l_stable = l_grid(k)
188	ENDIF
189	!
190	!-- Adjustment of the mixing length
191	IF ( wall_adjustment ) THEN
192	l(k,j) = MIN( wall_adjustment_factor * &
193	( zu(k) - zw(nzb_s_inner(j,i)) ), &
194	l_grid(k), l_stable )
195	ll(k,j) = MIN( wall_adjustment_factor * &
196	( zu(k) - zw(nzb_s_inner(j,i)) ), &
197	l_grid(k) )
198	ELSE
199	l(k,j) = MIN( l_grid(k), l_stable )
200	ll(k,j) = l_grid(k)
201	ENDIF
202
203	ENDDO
204	ENDDO
205
206	!
207	!-- Calculate the tendency terms
208	DO j = nys, nyn
209	DO k = nzb_s_inner(j,i)+1, nzt
210
211	dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) * &
212	e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
213
214	tend(k,j,i) = tend(k,j,i) &
215	+ ( &
216	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
217	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
218	) * ddx2 &
219	+ ( &
220	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
221	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
222	) * ddy2 &
223	+ ( &
224	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
225	* rho_air_zw(k) &
226	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
227	* rho_air_zw(k-1) &
228	) * ddzw(k) * drho_air(k) &
229	- dissipation(k,j)
230
231	ENDDO
232	ENDDO
233
234	!
235	!-- Store dissipation if needed for calculating the sgs particle
236	!-- velocities
237	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
238	collision_turbulence ) THEN
239	DO j = nys, nyn
240	DO k = nzb_s_inner(j,i)+1, nzt
241	diss(k,j,i) = dissipation(k,j)
242	ENDDO
243	ENDDO
244	ENDIF
245
246	ENDDO
247
248	ELSE
249
250	DO i = nxl, nxr
251	DO j = nys, nyn
252	DO k = nzb_s_inner(j,i)+1, nzt
253	!
254	!-- Calculate the mixing length (for dissipation)
255	dvar_dz = atmos_ocean_sign * &
256	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
257	IF ( dvar_dz > 0.0_wp ) THEN
258	l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
259	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
260	ELSE
261	l_stable = l_grid(k)
262	ENDIF
263	!
264	!-- Adjustment of the mixing length
265	IF ( wall_adjustment ) THEN
266	l(k,j) = MIN( wall_adjustment_factor * &
267	( zu(k) - zw(nzb_s_inner(j,i)) ), &
268	l_grid(k), l_stable )
269	ll(k,j) = MIN( wall_adjustment_factor * &
270	( zu(k) - zw(nzb_s_inner(j,i)) ), &
271	l_grid(k) )
272	ELSE
273	l(k,j) = MIN( l_grid(k), l_stable )
274	ll(k,j) = l_grid(k)
275	ENDIF
276
277	ENDDO
278	ENDDO
279
280	!
281	!-- Calculate the tendency terms
282	DO j = nys, nyn
283	DO k = nzb_s_inner(j,i)+1, nzt
284
285	dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) * &
286	e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
287
288	tend(k,j,i) = tend(k,j,i) &
289	+ ( &
290	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
291	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
292	) * ddx2 &
293	+ ( &
294	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
295	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
296	) * ddy2 &
297	+ ( &
298	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
299	* rho_air_zw(k) &
300	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
301	* rho_air_zw(k-1) &
302	) * ddzw(k) * drho_air(k) &
303	- dissipation(k,j)
304
305	ENDDO
306	ENDDO
307
308	!
309	!-- Store dissipation if needed for calculating the sgs particle
310	!-- velocities
311	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
312	collision_turbulence ) THEN
313	DO j = nys, nyn
314	DO k = nzb_s_inner(j,i)+1, nzt
315	diss(k,j,i) = dissipation(k,j)
316	ENDDO
317	ENDDO
318	ENDIF
319
320	ENDDO
321
322	ENDIF
323
324	!
325	!-- Boundary condition for dissipation
326	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
327	collision_turbulence ) THEN
328	DO i = nxl, nxr
329	DO j = nys, nyn
330	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
331	ENDDO
332	ENDDO
333	ENDIF
334
335	END SUBROUTINE diffusion_e
336
337
338	!------------------------------------------------------------------------------!
339	! Description:
340	! ------------
341	!> Call for all grid points - accelerator version
342	!------------------------------------------------------------------------------!
343	SUBROUTINE diffusion_e_acc( var, var_reference )
344
345	USE arrays_3d, &
346	ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw, &
347	drho_air, rho_air_zw
348
349	USE control_parameters, &
350	ONLY: atmos_ocean_sign, g, use_single_reference_value, &
351	wall_adjustment, wall_adjustment_factor
352
353	USE grid_variables, &
354	ONLY: ddx2, ddy2
355
356	USE indices, &
357	ONLY: i_left, i_right, j_north, j_south, nxlg, nxrg, nyng, nysg, &
358	nzb, nzb_s_inner, nzt
359
360	USE kinds
361
362	USE microphysics_mod, &
363	ONLY: collision_turbulence
364
365	USE particle_attributes, &
366	ONLY: use_sgs_for_particles, wang_kernel
367
368	IMPLICIT NONE
369
370	INTEGER(iwp) :: i !<
371	INTEGER(iwp) :: j !<
372	INTEGER(iwp) :: k !<
373	REAL(wp) :: dissipation !<
374	REAL(wp) :: dvar_dz !<
375	REAL(wp) :: l !<
376	REAL(wp) :: ll !<
377	REAL(wp) :: l_stable !<
378	REAL(wp) :: var_reference !<
379
380	#if defined( __nopointer )
381	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !<
382	#else
383	REAL(wp), DIMENSION(:,:,:), POINTER :: var !<
384	#endif
385
386
387	!
388	!-- This if clause must be outside the k-loop because otherwise
389	!-- runtime errors occur with -C hopt on NEC
390	IF ( use_single_reference_value ) THEN
391
392	!$acc kernels present( ddzu, ddzw, dd2zu, diss, e, km, l_grid ) &
393	!$acc present( nzb_s_inner, tend, var, zu, zw )
394	DO i = i_left, i_right
395	DO j = j_south, j_north
396	DO k = 1, nzt
397
398	IF ( k > nzb_s_inner(j,i) ) THEN
399	!
400	!-- Calculate the mixing length (for dissipation)
401	dvar_dz = atmos_ocean_sign * &
402	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
403	IF ( dvar_dz > 0.0_wp ) THEN
404	l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
405	SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
406	ELSE
407	l_stable = l_grid(k)
408	ENDIF
409	!
410	!-- Adjustment of the mixing length
411	IF ( wall_adjustment ) THEN
412	l = MIN( wall_adjustment_factor * &
413	( zu(k) - zw(nzb_s_inner(j,i)) ), &
414	l_grid(k), l_stable )
415	ll = MIN( wall_adjustment_factor * &
416	( zu(k) - zw(nzb_s_inner(j,i)) ), &
417	l_grid(k) )
418	ELSE
419	l = MIN( l_grid(k), l_stable )
420	ll = l_grid(k)
421	ENDIF
422	!
423	!-- Calculate the tendency terms
424	dissipation = ( 0.19_wp + 0.74_wp * l / ll ) * &
425	e(k,j,i) * SQRT( e(k,j,i) ) / l
426
427	tend(k,j,i) = tend(k,j,i) &
428	+ ( &
429	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
430	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
431	) * ddx2 &
432	+ ( &
433	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
434	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
435	) * ddy2 &
436	+ ( &
437	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
438	* rho_air_zw(k) &
439	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
440	* rho_air_zw(k-1) &
441	) * ddzw(k) * drho_air(k) &
442	- dissipation
443
444	!
445	!-- Store dissipation if needed for calculating the sgs particle
446	!-- velocities
447	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
448	collision_turbulence ) THEN
449	diss(k,j,i) = dissipation
450	ENDIF
451
452	ENDIF
453
454	ENDDO
455	ENDDO
456	ENDDO
457	!$acc end kernels
458
459	ELSE
460
461	!$acc kernels present( ddzu, ddzw, dd2zu, diss, e, km, l_grid ) &
462	!$acc present( nzb_s_inner, tend, var, zu, zw )
463	DO i = i_left, i_right
464	DO j = j_south, j_north
465	DO k = 1, nzt
466
467	IF ( k > nzb_s_inner(j,i) ) THEN
468	!
469	!-- Calculate the mixing length (for dissipation)
470	dvar_dz = atmos_ocean_sign * &
471	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
472	IF ( dvar_dz > 0.0_wp ) THEN
473	l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
474	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
475	ELSE
476	l_stable = l_grid(k)
477	ENDIF
478	!
479	!-- Adjustment of the mixing length
480	IF ( wall_adjustment ) THEN
481	l = MIN( wall_adjustment_factor * &
482	( zu(k) - zw(nzb_s_inner(j,i)) ), &
483	l_grid(k), l_stable )
484	ll = MIN( wall_adjustment_factor * &
485	( zu(k) - zw(nzb_s_inner(j,i)) ), &
486	l_grid(k) )
487	ELSE
488	l = MIN( l_grid(k), l_stable )
489	ll = l_grid(k)
490	ENDIF
491	!
492	!-- Calculate the tendency terms
493	dissipation = ( 0.19_wp + 0.74_wp * l / ll ) * &
494	e(k,j,i) * SQRT( e(k,j,i) ) / l
495
496	tend(k,j,i) = tend(k,j,i) &
497	+ ( &
498	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
499	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
500	) * ddx2 &
501	+ ( &
502	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
503	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
504	) * ddy2 &
505	+ ( &
506	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
507	* rho_air_zw(k) &
508	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
509	* rho_air_zw(k-1) &
510	) * ddzw(k) * drho_air(k) &
511	- dissipation
512
513	!
514	!-- Store dissipation if needed for calculating the sgs
515	!-- particle velocities
516	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
517	collision_turbulence ) THEN
518	diss(k,j,i) = dissipation
519	ENDIF
520
521	ENDIF
522
523	ENDDO
524	ENDDO
525	ENDDO
526	!$acc end kernels
527
528	ENDIF
529
530	!
531	!-- Boundary condition for dissipation
532	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
533	collision_turbulence ) THEN
534	!$acc kernels present( diss, nzb_s_inner )
535	DO i = i_left, i_right
536	DO j = j_south, j_north
537	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
538	ENDDO
539	ENDDO
540	!$acc end kernels
541	ENDIF
542
543	END SUBROUTINE diffusion_e_acc
544
545
546	!------------------------------------------------------------------------------!
547	! Description:
548	! ------------
549	!> Call for grid point i,j
550	!------------------------------------------------------------------------------!
551	SUBROUTINE diffusion_e_ij( i, j, var, var_reference )
552
553	USE arrays_3d, &
554	ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw, &
555	drho_air, rho_air_zw
556
557	USE control_parameters, &
558	ONLY: atmos_ocean_sign, g, use_single_reference_value, &
559	wall_adjustment, wall_adjustment_factor
560
561	USE grid_variables, &
562	ONLY: ddx2, ddy2
563
564	USE indices, &
565	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzb_s_inner, nzt
566
567	USE kinds
568
569	USE microphysics_mod, &
570	ONLY: collision_turbulence
571
572	USE particle_attributes, &
573	ONLY: use_sgs_for_particles, wang_kernel
574
575	IMPLICIT NONE
576
577	INTEGER(iwp) :: i !<
578	INTEGER(iwp) :: j !<
579	INTEGER(iwp) :: k !<
580	REAL(wp) :: dvar_dz !<
581	REAL(wp) :: l_stable !<
582	REAL(wp) :: var_reference !<
583
584	#if defined( __nopointer )
585	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !<
586	#else
587	REAL(wp), DIMENSION(:,:,:), POINTER :: var !<
588	#endif
589	REAL(wp), DIMENSION(nzb+1:nzt) :: dissipation !<
590	REAL(wp), DIMENSION(nzb+1:nzt) :: l !<
591	REAL(wp), DIMENSION(nzb+1:nzt) :: ll !<
592
593
594	!
595	!-- Calculate the mixing length (for dissipation)
596	DO k = nzb_s_inner(j,i)+1, nzt
597	dvar_dz = atmos_ocean_sign * &
598	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
599	IF ( dvar_dz > 0.0_wp ) THEN
600	IF ( use_single_reference_value ) THEN
601	l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
602	SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
603	ELSE
604	l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
605	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
606	ENDIF
607	ELSE
608	l_stable = l_grid(k)
609	ENDIF
610	!
611	!-- Adjustment of the mixing length
612	IF ( wall_adjustment ) THEN
613	l(k) = MIN( wall_adjustment_factor * &
614	( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k), &
615	l_stable )
616	ll(k) = MIN( wall_adjustment_factor * &
617	( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k) )
618	ELSE
619	l(k) = MIN( l_grid(k), l_stable )
620	ll(k) = l_grid(k)
621	ENDIF
622	!
623	!-- Calculate the tendency term
624	dissipation(k) = ( 0.19_wp + 0.74_wp * l(k) / ll(k) ) * e(k,j,i) * &
625	SQRT( e(k,j,i) ) / l(k)
626
627	tend(k,j,i) = tend(k,j,i) &
628	+ ( &
629	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
630	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
631	) * ddx2 &
632	+ ( &
633	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
634	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
635	) * ddy2 &
636	+ ( &
637	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
638	* rho_air_zw(k) &
639	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
640	* rho_air_zw(k-1) &
641	) * ddzw(k) * drho_air(k) &
642	- dissipation(k)
643
644	ENDDO
645
646	!
647	!-- Store dissipation if needed for calculating the sgs particle velocities
648	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
649	collision_turbulence ) THEN
650	DO k = nzb_s_inner(j,i)+1, nzt
651	diss(k,j,i) = dissipation(k)
652	ENDDO
653	!
654	!-- Boundary condition for dissipation
655	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
656	ENDIF
657
658	END SUBROUTINE diffusion_e_ij
659
660	END MODULE diffusion_e_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |