Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |