Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |