Home

Context Navigation

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

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