Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |