Home

Context Navigation

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

Last change on this file since 1095 was 1066, checked in by hoffmann, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 21.5 KB

Line
1	MODULE diffusion_e_mod
2
3	!--------------------------------------------------------------------------------!
4	! This file is part of PALM.
5	!
6	! PALM is free software: you can redistribute it and/or modify it under the terms
7	! of the GNU General Public License as published by the Free Software Foundation,
8	! either version 3 of the License, or (at your option) any later 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-2012 Leibniz University Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	! Former revisions:
24	! -----------------
25	! $Id: diffusion_e.f90 1066 2012-11-22 17:52:43Z raasch $
26	!
27	! 1065 2012-11-22 17:42:36Z hoffmann
28	! Enabled the claculation of diss in case of turbulence = .TRUE. (parameterized
29	! effects of turbulence on autoconversion and accretion in two-moments cloud
30	! physics scheme).
31	!
32	! 1036 2012-10-22 13:43:42Z raasch
33	! code put under GPL (PALM 3.9)
34	!
35	! 1015 2012-09-27 09:23:24Z raasch
36	! accelerator version (*_acc) added,
37	! adjustment of mixing length to the Prandtl mixing length at first grid point
38	! above ground removed
39	!
40	! 1010 2012-09-20 07:59:54Z raasch
41	! cpp switch __nopointer added for pointer free version
42	!
43	! 1001 2012-09-13 14:08:46Z raasch
44	! most arrays comunicated by module instead of parameter list
45	!
46	! 825 2012-02-19 03:03:44Z raasch
47	! wang_collision_kernel renamed wang_kernel
48	!
49	! 790 2011-11-29 03:11:20Z raasch
50	! diss is also calculated in case that the Wang kernel is used
51	!
52	! 667 2010-12-23 12:06:00Z suehring/gryschka
53	! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
54	!
55	! 97 2007-06-21 08:23:15Z raasch
56	! Adjustment of mixing length calculation for the ocean version. zw added to
57	! argument list.
58	! This is also a bugfix, because the height above the topography is now
59	! used instead of the height above level k=0.
60	! theta renamed var, dpt_dz renamed dvar_dz, +new argument var_reference
61	! use_pt_reference renamed use_reference
62	!
63	! 65 2007-03-13 12:11:43Z raasch
64	! Reference temperature pt_reference can be used in buoyancy term
65	!
66	! 20 2007-02-26 00:12:32Z raasch
67	! Bugfix: ddzw dimensioned 1:nzt"+1"
68	! Calculation extended for gridpoint nzt
69	!
70	! RCS Log replace by Id keyword, revision history cleaned up
71	!
72	! Revision 1.18 2006/08/04 14:29:43 raasch
73	! dissipation is stored in extra array diss if needed later on for calculating
74	! the sgs particle velocities
75	!
76	! Revision 1.1 1997/09/19 07:40:24 raasch
77	! Initial revision
78	!
79	!
80	! Description:
81	! ------------
82	! Diffusion- and dissipation terms for the TKE
83	!------------------------------------------------------------------------------!
84
85	PRIVATE
86	PUBLIC diffusion_e, diffusion_e_acc
87
88
89	INTERFACE diffusion_e
90	MODULE PROCEDURE diffusion_e
91	MODULE PROCEDURE diffusion_e_ij
92	END INTERFACE diffusion_e
93
94	INTERFACE diffusion_e_acc
95	MODULE PROCEDURE diffusion_e_acc
96	END INTERFACE diffusion_e_acc
97
98	CONTAINS
99
100
101	!------------------------------------------------------------------------------!
102	! Call for all grid points
103	!------------------------------------------------------------------------------!
104	SUBROUTINE diffusion_e( var, var_reference )
105
106	USE arrays_3d
107	USE control_parameters
108	USE grid_variables
109	USE indices
110	USE particle_attributes
111
112	IMPLICIT NONE
113
114	INTEGER :: i, j, k
115	REAL :: dvar_dz, l_stable, var_reference
116
117	#if defined( __nopointer )
118	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var
119	#else
120	REAL, DIMENSION(:,:,:), POINTER :: var
121	#endif
122	REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: dissipation, l, ll
123
124
125	!
126	!-- This if clause must be outside the k-loop because otherwise
127	!-- runtime errors occur with -C hopt on NEC
128	IF ( use_reference ) THEN
129
130	DO i = nxl, nxr
131	DO j = nys, nyn
132	DO k = nzb_s_inner(j,i)+1, nzt
133	!
134	!-- Calculate the mixing length (for dissipation)
135	dvar_dz = atmos_ocean_sign * &
136	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
137	IF ( dvar_dz > 0.0 ) THEN
138	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
139	SQRT( g / var_reference * dvar_dz ) + 1E-5
140	ELSE
141	l_stable = l_grid(k)
142	ENDIF
143	!
144	!-- Adjustment of the mixing length
145	IF ( wall_adjustment ) THEN
146	l(k,j) = MIN( wall_adjustment_factor * &
147	( zu(k) - zw(nzb_s_inner(j,i)) ), &
148	l_grid(k), l_stable )
149	ll(k,j) = MIN( wall_adjustment_factor * &
150	( zu(k) - zw(nzb_s_inner(j,i)) ), &
151	l_grid(k) )
152	ELSE
153	l(k,j) = MIN( l_grid(k), l_stable )
154	ll(k,j) = l_grid(k)
155	ENDIF
156
157	ENDDO
158	ENDDO
159
160	!
161	!-- Calculate the tendency terms
162	DO j = nys, nyn
163	DO k = nzb_s_inner(j,i)+1, nzt
164
165	dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
166	e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
167
168	tend(k,j,i) = tend(k,j,i) &
169	+ ( &
170	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
171	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
172	) * ddx2 &
173	+ ( &
174	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
175	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
176	) * ddy2 &
177	+ ( &
178	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
179	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
180	) * ddzw(k) &
181	- dissipation(k,j)
182
183	ENDDO
184	ENDDO
185
186	!
187	!-- Store dissipation if needed for calculating the sgs particle
188	!-- velocities
189	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
190	turbulence ) THEN
191	DO j = nys, nyn
192	DO k = nzb_s_inner(j,i)+1, nzt
193	diss(k,j,i) = dissipation(k,j)
194	ENDDO
195	ENDDO
196	ENDIF
197
198	ENDDO
199
200	ELSE
201
202	DO i = nxl, nxr
203	DO j = nys, nyn
204	DO k = nzb_s_inner(j,i)+1, nzt
205	!
206	!-- Calculate the mixing length (for dissipation)
207	dvar_dz = atmos_ocean_sign * &
208	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
209	IF ( dvar_dz > 0.0 ) THEN
210	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
211	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
212	ELSE
213	l_stable = l_grid(k)
214	ENDIF
215	!
216	!-- Adjustment of the mixing length
217	IF ( wall_adjustment ) THEN
218	l(k,j) = MIN( wall_adjustment_factor * &
219	( zu(k) - zw(nzb_s_inner(j,i)) ), &
220	l_grid(k), l_stable )
221	ll(k,j) = MIN( wall_adjustment_factor * &
222	( zu(k) - zw(nzb_s_inner(j,i)) ), &
223	l_grid(k) )
224	ELSE
225	l(k,j) = MIN( l_grid(k), l_stable )
226	ll(k,j) = l_grid(k)
227	ENDIF
228
229	ENDDO
230	ENDDO
231
232	!
233	!-- Calculate the tendency terms
234	DO j = nys, nyn
235	DO k = nzb_s_inner(j,i)+1, nzt
236
237	dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
238	e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
239
240	tend(k,j,i) = tend(k,j,i) &
241	+ ( &
242	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
243	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
244	) * ddx2 &
245	+ ( &
246	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
247	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
248	) * ddy2 &
249	+ ( &
250	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
251	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
252	) * ddzw(k) &
253	- dissipation(k,j)
254
255	ENDDO
256	ENDDO
257
258	!
259	!-- Store dissipation if needed for calculating the sgs particle
260	!-- velocities
261	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
262	turbulence ) THEN
263	DO j = nys, nyn
264	DO k = nzb_s_inner(j,i)+1, nzt
265	diss(k,j,i) = dissipation(k,j)
266	ENDDO
267	ENDDO
268	ENDIF
269
270	ENDDO
271
272	ENDIF
273
274	!
275	!-- Boundary condition for dissipation
276	IF ( use_sgs_for_particles .OR. wang_kernel .OR. turbulence ) THEN
277	DO i = nxl, nxr
278	DO j = nys, nyn
279	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
280	ENDDO
281	ENDDO
282	ENDIF
283
284	END SUBROUTINE diffusion_e
285
286
287	!------------------------------------------------------------------------------!
288	! Call for all grid points - accelerator version
289	!------------------------------------------------------------------------------!
290	SUBROUTINE diffusion_e_acc( var, var_reference )
291
292	USE arrays_3d
293	USE control_parameters
294	USE grid_variables
295	USE indices
296	USE particle_attributes
297
298	IMPLICIT NONE
299
300	INTEGER :: i, j, k
301	REAL :: dissipation, dvar_dz, l, ll, l_stable, var_reference
302
303	#if defined( __nopointer )
304	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var
305	#else
306	REAL, DIMENSION(:,:,:), POINTER :: var
307	#endif
308
309
310	!
311	!-- This if clause must be outside the k-loop because otherwise
312	!-- runtime errors occur with -C hopt on NEC
313	IF ( use_reference ) THEN
314	STOP '+++ use_reference in diffusion_e not implemented'
315	! DO i = nxl, nxr
316	! DO j = nys, nyn
317	! DO k = nzb_s_inner(j,i)+1, nzt
318	!
319	!-- Calculate the mixing length (for dissipation)
320	! dvar_dz = atmos_ocean_sign * &
321	! ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
322	! IF ( dvar_dz > 0.0 ) THEN
323	! l_stable = 0.76 * SQRT( e(k,j,i) ) / &
324	! SQRT( g / var_reference * dvar_dz ) + 1E-5
325	! ELSE
326	! l_stable = l_grid(k)
327	! ENDIF
328	!
329	!-- Adjustment of the mixing length
330	! IF ( wall_adjustment ) THEN
331	! l(k,j) = MIN( wall_adjustment_factor * &
332	! ( zu(k) - zw(nzb_s_inner(j,i)) ), &
333	! l_grid(k), l_stable )
334	! ll(k,j) = MIN( wall_adjustment_factor * &
335	! ( zu(k) - zw(nzb_s_inner(j,i)) ), &
336	! l_grid(k) )
337	! ELSE
338	! l(k,j) = MIN( l_grid(k), l_stable )
339	! ll(k,j) = l_grid(k)
340	! ENDIF
341	!
342	! ENDDO
343	! ENDDO
344	!
345	!
346	!-- Calculate the tendency terms
347	! DO j = nys, nyn
348	! DO k = nzb_s_inner(j,i)+1, nzt
349	!
350	! dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
351	! e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
352	!
353	! tend(k,j,i) = tend(k,j,i) &
354	! + ( &
355	! ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
356	! - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
357	! ) * ddx2 &
358	! + ( &
359	! ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
360	! - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
361	! ) * ddy2 &
362	! + ( &
363	! ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
364	! - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
365	! ) * ddzw(k) &
366	! - dissipation(k,j)
367	!
368	! ENDDO
369	! ENDDO
370	!
371	!
372	!-- Store dissipation if needed for calculating the sgs particle
373	!-- velocities
374	! IF ( use_sgs_for_particles .OR. wang_kernel ) THEN
375	! DO j = nys, nyn
376	! DO k = nzb_s_inner(j,i)+1, nzt
377	! diss(k,j,i) = dissipation(k,j)
378	! ENDDO
379	! ENDDO
380	! ENDIF
381	!
382	! ENDDO
383	!
384	ELSE
385
386	!$acc kernels present( ddzu, ddzw, dd2zu, diss, e, km, l_grid ) &
387	!$acc present( nzb_s_inner, rif, tend, var, zu, zw )
388	!$acc loop
389	DO i = nxl, nxr
390	DO j = nys, nyn
391	!$acc loop vector( 32 )
392	DO k = 1, nzt
393
394	IF ( k > nzb_s_inner(j,i) ) THEN
395	!
396	!-- Calculate the mixing length (for dissipation)
397	dvar_dz = atmos_ocean_sign * &
398	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
399	IF ( dvar_dz > 0.0 ) THEN
400	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
401	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
402	ELSE
403	l_stable = l_grid(k)
404	ENDIF
405	!
406	!-- Adjustment of the mixing length
407	IF ( wall_adjustment ) THEN
408	l = MIN( wall_adjustment_factor * &
409	( zu(k) - zw(nzb_s_inner(j,i)) ), &
410	l_grid(k), l_stable )
411	ll = MIN( wall_adjustment_factor * &
412	( zu(k) - zw(nzb_s_inner(j,i)) ), &
413	l_grid(k) )
414	ELSE
415	l = MIN( l_grid(k), l_stable )
416	ll = l_grid(k)
417	ENDIF
418	!
419	!-- Calculate the tendency terms
420	dissipation = ( 0.19 + 0.74 * l / ll ) * &
421	e(k,j,i) * SQRT( e(k,j,i) ) / l
422
423	tend(k,j,i) = tend(k,j,i) &
424	+ ( &
425	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
426	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
427	) * ddx2 &
428	+ ( &
429	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
430	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
431	) * ddy2 &
432	+ ( &
433	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
434	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
435	) * ddzw(k) &
436	- dissipation
437
438	!
439	!-- Store dissipation if needed for calculating the sgs
440	!-- particle velocities
441	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
442	turbulence ) THEN
443	diss(k,j,i) = dissipation
444	ENDIF
445
446	ENDIF
447
448	ENDDO
449	ENDDO
450	ENDDO
451	!$acc end kernels
452
453	ENDIF
454
455	!
456	!-- Boundary condition for dissipation
457	IF ( use_sgs_for_particles .OR. wang_kernel .OR. turbulence ) THEN
458	!$acc kernels present( diss, nzb_s_inner )
459	!$acc loop
460	DO i = nxl, nxr
461	!$acc loop vector( 32 )
462	DO j = nys, nyn
463	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
464	ENDDO
465	ENDDO
466	!$acc end kernels
467	ENDIF
468
469	END SUBROUTINE diffusion_e_acc
470
471
472	!------------------------------------------------------------------------------!
473	! Call for grid point i,j
474	!------------------------------------------------------------------------------!
475	SUBROUTINE diffusion_e_ij( i, j, var, var_reference )
476
477	USE arrays_3d
478	USE control_parameters
479	USE grid_variables
480	USE indices
481	USE particle_attributes
482
483	IMPLICIT NONE
484
485	INTEGER :: i, j, k
486	REAL :: dvar_dz, l_stable, var_reference
487
488	#if defined( __nopointer )
489	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var
490	#else
491	REAL, DIMENSION(:,:,:), POINTER :: var
492	#endif
493	REAL, DIMENSION(nzb+1:nzt) :: dissipation, l, ll
494
495
496	!
497	!-- Calculate the mixing length (for dissipation)
498	DO k = nzb_s_inner(j,i)+1, nzt
499	dvar_dz = atmos_ocean_sign * &
500	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
501	IF ( dvar_dz > 0.0 ) THEN
502	IF ( use_reference ) THEN
503	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
504	SQRT( g / var_reference * dvar_dz ) + 1E-5
505	ELSE
506	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
507	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
508	ENDIF
509	ELSE
510	l_stable = l_grid(k)
511	ENDIF
512	!
513	!-- Adjustment of the mixing length
514	IF ( wall_adjustment ) THEN
515	l(k) = MIN( wall_adjustment_factor * &
516	( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k), &
517	l_stable )
518	ll(k) = MIN( wall_adjustment_factor * &
519	( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k) )
520	ELSE
521	l(k) = MIN( l_grid(k), l_stable )
522	ll(k) = l_grid(k)
523	ENDIF
524	!
525	!-- Calculate the tendency term
526	dissipation(k) = ( 0.19 + 0.74 * l(k) / ll(k) ) * e(k,j,i) * &
527	SQRT( e(k,j,i) ) / l(k)
528
529	tend(k,j,i) = tend(k,j,i) &
530	+ ( &
531	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
532	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
533	) * ddx2 &
534	+ ( &
535	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
536	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
537	) * ddy2 &
538	+ ( &
539	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
540	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
541	) * ddzw(k) &
542	- dissipation(k)
543
544	ENDDO
545
546	!
547	!-- Store dissipation if needed for calculating the sgs particle velocities
548	IF ( use_sgs_for_particles .OR. wang_kernel .OR. turbulence ) THEN
549	DO k = nzb_s_inner(j,i)+1, nzt
550	diss(k,j,i) = dissipation(k)
551	ENDDO
552	!
553	!-- Boundary condition for dissipation
554	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
555	ENDIF
556
557	END SUBROUTINE diffusion_e_ij
558
559	END MODULE diffusion_e_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |