Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |