Home

Context Navigation

diffusion_e.f90 @ 1171

Last change on this file since 1171 was 1171, checked in by raasch, 11 years ago

New:
---

use_reference-case activated in accelerator version. (buoyancy, diffusion_e)
new option -e which defines the execution command to be used to run PALM,
compiler options for pgi/openacc added (palm_simple_run)
parameter sets for openACC benchmarks added (trunk/EXAMPLES/benchmark_acc)

Changed:

split of prognostic_equations deactivated (time_integration)

Errors:

bugfix: diss array is allocated with full size if accelerator boards are used (init_3d_model)

Property svn:keywords set to Id

File size: 21.7 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	! use_reference-case activated in accelerator version
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: diffusion_e.f90 1171 2013-05-30 11:27:45Z raasch $
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
320	!$acc kernels present( ddzu, ddzw, dd2zu, diss, e, km, l_grid ) &
321	!$acc present( nzb_s_inner, rif, tend, var, zu, zw )
322	!$acc loop
323	DO i = i_left, i_right
324	DO j = j_south, j_north
325	!$acc loop vector( 32 )
326	DO k = 1, nzt
327
328	IF ( k > nzb_s_inner(j,i) ) THEN
329	!
330	!-- Calculate the mixing length (for dissipation)
331	dvar_dz = atmos_ocean_sign * &
332	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
333	IF ( dvar_dz > 0.0 ) THEN
334	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
335	SQRT( g / var_reference * dvar_dz ) + 1E-5
336	ELSE
337	l_stable = l_grid(k)
338	ENDIF
339	!
340	!-- Adjustment of the mixing length
341	IF ( wall_adjustment ) THEN
342	l = MIN( wall_adjustment_factor * &
343	( zu(k) - zw(nzb_s_inner(j,i)) ), &
344	l_grid(k), l_stable )
345	ll = MIN( wall_adjustment_factor * &
346	( zu(k) - zw(nzb_s_inner(j,i)) ), &
347	l_grid(k) )
348	ELSE
349	l = MIN( l_grid(k), l_stable )
350	ll = l_grid(k)
351	ENDIF
352	!
353	!-- Calculate the tendency terms
354	dissipation = ( 0.19 + 0.74 * l / ll ) * &
355	e(k,j,i) * SQRT( e(k,j,i) ) / l
356
357	tend(k,j,i) = tend(k,j,i) &
358	+ ( &
359	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
360	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
361	) * ddx2 &
362	+ ( &
363	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
364	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
365	) * ddy2 &
366	+ ( &
367	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
368	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
369	) * ddzw(k) &
370	- dissipation
371
372	!
373	!-- Store dissipation if needed for calculating the sgs particle
374	!-- velocities
375	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
376	turbulence ) THEN
377	diss(k,j,i) = dissipation
378	ENDIF
379
380	ENDIF
381
382	ENDDO
383	ENDDO
384	ENDDO
385	!$acc end kernels
386
387	ELSE
388
389	!$acc kernels present( ddzu, ddzw, dd2zu, diss, e, km, l_grid ) &
390	!$acc present( nzb_s_inner, rif, tend, var, zu, zw )
391	!$acc loop
392	DO i = i_left, i_right
393	DO j = j_south, j_north
394	!$acc loop vector( 32 )
395	DO k = 1, nzt
396
397	IF ( k > nzb_s_inner(j,i) ) THEN
398	!
399	!-- Calculate the mixing length (for dissipation)
400	dvar_dz = atmos_ocean_sign * &
401	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
402	IF ( dvar_dz > 0.0 ) THEN
403	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
404	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
405	ELSE
406	l_stable = l_grid(k)
407	ENDIF
408	!
409	!-- Adjustment of the mixing length
410	IF ( wall_adjustment ) THEN
411	l = MIN( wall_adjustment_factor * &
412	( zu(k) - zw(nzb_s_inner(j,i)) ), &
413	l_grid(k), l_stable )
414	ll = MIN( wall_adjustment_factor * &
415	( zu(k) - zw(nzb_s_inner(j,i)) ), &
416	l_grid(k) )
417	ELSE
418	l = MIN( l_grid(k), l_stable )
419	ll = l_grid(k)
420	ENDIF
421	!
422	!-- Calculate the tendency terms
423	dissipation = ( 0.19 + 0.74 * l / ll ) * &
424	e(k,j,i) * SQRT( e(k,j,i) ) / l
425
426	tend(k,j,i) = tend(k,j,i) &
427	+ ( &
428	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
429	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
430	) * ddx2 &
431	+ ( &
432	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
433	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
434	) * ddy2 &
435	+ ( &
436	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
437	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
438	) * ddzw(k) &
439	- dissipation
440
441	!
442	!-- Store dissipation if needed for calculating the sgs
443	!-- particle velocities
444	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
445	turbulence ) THEN
446	diss(k,j,i) = dissipation
447	ENDIF
448
449	ENDIF
450
451	ENDDO
452	ENDDO
453	ENDDO
454	!$acc end kernels
455
456	ENDIF
457
458	!
459	!-- Boundary condition for dissipation
460	IF ( use_sgs_for_particles .OR. wang_kernel .OR. turbulence ) THEN
461	!$acc kernels present( diss, nzb_s_inner )
462	!$acc loop
463	DO i = i_left, i_right
464	!$acc loop vector( 32 )
465	DO j = j_south, j_north
466	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
467	ENDDO
468	ENDDO
469	!$acc end kernels
470	ENDIF
471
472	END SUBROUTINE diffusion_e_acc
473
474
475	!------------------------------------------------------------------------------!
476	! Call for grid point i,j
477	!------------------------------------------------------------------------------!
478	SUBROUTINE diffusion_e_ij( i, j, var, var_reference )
479
480	USE arrays_3d
481	USE control_parameters
482	USE grid_variables
483	USE indices
484	USE particle_attributes
485
486	IMPLICIT NONE
487
488	INTEGER :: i, j, k
489	REAL :: dvar_dz, l_stable, var_reference
490
491	#if defined( __nopointer )
492	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var
493	#else
494	REAL, DIMENSION(:,:,:), POINTER :: var
495	#endif
496	REAL, DIMENSION(nzb+1:nzt) :: dissipation, l, ll
497
498
499	!
500	!-- Calculate the mixing length (for dissipation)
501	DO k = nzb_s_inner(j,i)+1, nzt
502	dvar_dz = atmos_ocean_sign * &
503	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
504	IF ( dvar_dz > 0.0 ) THEN
505	IF ( use_reference ) THEN
506	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
507	SQRT( g / var_reference * dvar_dz ) + 1E-5
508	ELSE
509	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
510	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
511	ENDIF
512	ELSE
513	l_stable = l_grid(k)
514	ENDIF
515	!
516	!-- Adjustment of the mixing length
517	IF ( wall_adjustment ) THEN
518	l(k) = MIN( wall_adjustment_factor * &
519	( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k), &
520	l_stable )
521	ll(k) = MIN( wall_adjustment_factor * &
522	( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k) )
523	ELSE
524	l(k) = MIN( l_grid(k), l_stable )
525	ll(k) = l_grid(k)
526	ENDIF
527	!
528	!-- Calculate the tendency term
529	dissipation(k) = ( 0.19 + 0.74 * l(k) / ll(k) ) * e(k,j,i) * &
530	SQRT( e(k,j,i) ) / l(k)
531
532	tend(k,j,i) = tend(k,j,i) &
533	+ ( &
534	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
535	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
536	) * ddx2 &
537	+ ( &
538	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
539	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
540	) * ddy2 &
541	+ ( &
542	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
543	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
544	) * ddzw(k) &
545	- dissipation(k)
546
547	ENDDO
548
549	!
550	!-- Store dissipation if needed for calculating the sgs particle velocities
551	IF ( use_sgs_for_particles .OR. wang_kernel .OR. turbulence ) THEN
552	DO k = nzb_s_inner(j,i)+1, nzt
553	diss(k,j,i) = dissipation(k)
554	ENDDO
555	!
556	!-- Boundary condition for dissipation
557	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
558	ENDIF
559
560	END SUBROUTINE diffusion_e_ij
561
562	END MODULE diffusion_e_mod

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

Context Navigation

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

Download in other formats: