Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |