Home

Context Navigation

diffusion_e.f90 @ 1257

Last change on this file since 1257 was 1257, checked in by raasch, 10 years ago

New:
---

openACC porting of timestep calculation
(modules, timestep, time_integration)

Changed:

openACC loop directives and vector clauses removed (because they do not give any performance improvement with PGI
compiler versions > 13.6)
(advec_ws, buoyancy, coriolis, diffusion_e, diffusion_s, diffusion_u, diffusion_v, diffusion_w, diffusivities, exchange_horiz, fft_xy, pres, production_e, transpose, tridia_solver, wall_fluxes)

openACC loop independent clauses added
(boundary_conds, prandtl_fluxes, pres)

openACC declare create statements moved after FORTRAN declaration statement
(diffusion_u, diffusion_v, diffusion_w, fft_xy, poisfft, production_e, tridia_solver)

openACC end parallel replaced by end parallel loop
(flow_statistics, pres)

openACC "kernels do" replaced by "kernels loop"
(prandtl_fluxes)

output format for theta* changed to avoid output of *
(run_control)

Errors:

bugfix for calculation of advective timestep (old version may cause wrong timesteps in case of
vertixcally stretched grids)
Attention: standard run-control output has changed!
(timestep)

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-2012 Leibniz University Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	! openacc loop and loop vector clauses removed
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: diffusion_e.f90 1257 2013-11-08 15:18:40Z 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	DO i = i_left, i_right
329	DO j = j_south, j_north
330	DO k = 1, nzt
331
332	IF ( k > nzb_s_inner(j,i) ) THEN
333	!
334	!-- Calculate the mixing length (for dissipation)
335	dvar_dz = atmos_ocean_sign * &
336	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
337	IF ( dvar_dz > 0.0 ) THEN
338	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
339	SQRT( g / var_reference * dvar_dz ) + 1E-5
340	ELSE
341	l_stable = l_grid(k)
342	ENDIF
343	!
344	!-- Adjustment of the mixing length
345	IF ( wall_adjustment ) THEN
346	l = MIN( wall_adjustment_factor * &
347	( zu(k) - zw(nzb_s_inner(j,i)) ), &
348	l_grid(k), l_stable )
349	ll = MIN( wall_adjustment_factor * &
350	( zu(k) - zw(nzb_s_inner(j,i)) ), &
351	l_grid(k) )
352	ELSE
353	l = MIN( l_grid(k), l_stable )
354	ll = l_grid(k)
355	ENDIF
356	!
357	!-- Calculate the tendency terms
358	dissipation = ( 0.19 + 0.74 * l / ll ) * &
359	e(k,j,i) * SQRT( e(k,j,i) ) / l
360
361	tend(k,j,i) = tend(k,j,i) &
362	+ ( &
363	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
364	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
365	) * ddx2 &
366	+ ( &
367	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
368	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
369	) * ddy2 &
370	+ ( &
371	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
372	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
373	) * ddzw(k) &
374	- dissipation
375
376	!
377	!-- Store dissipation if needed for calculating the sgs particle
378	!-- velocities
379	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
380	turbulence ) THEN
381	diss(k,j,i) = dissipation
382	ENDIF
383
384	ENDIF
385
386	ENDDO
387	ENDDO
388	ENDDO
389	!$acc end kernels
390
391	ELSE
392
393	!$acc kernels present( ddzu, ddzw, dd2zu, diss, e, km, l_grid ) &
394	!$acc present( nzb_s_inner, rif, tend, var, zu, zw )
395	DO i = i_left, i_right
396	DO j = j_south, j_north
397	DO k = 1, nzt
398
399	IF ( k > nzb_s_inner(j,i) ) THEN
400	!
401	!-- Calculate the mixing length (for dissipation)
402	dvar_dz = atmos_ocean_sign * &
403	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
404	IF ( dvar_dz > 0.0 ) THEN
405	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
406	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
407	ELSE
408	l_stable = l_grid(k)
409	ENDIF
410	!
411	!-- Adjustment of the mixing length
412	IF ( wall_adjustment ) THEN
413	l = MIN( wall_adjustment_factor * &
414	( zu(k) - zw(nzb_s_inner(j,i)) ), &
415	l_grid(k), l_stable )
416	ll = MIN( wall_adjustment_factor * &
417	( zu(k) - zw(nzb_s_inner(j,i)) ), &
418	l_grid(k) )
419	ELSE
420	l = MIN( l_grid(k), l_stable )
421	ll = l_grid(k)
422	ENDIF
423	!
424	!-- Calculate the tendency terms
425	dissipation = ( 0.19 + 0.74 * l / ll ) * &
426	e(k,j,i) * SQRT( e(k,j,i) ) / l
427
428	tend(k,j,i) = tend(k,j,i) &
429	+ ( &
430	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
431	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
432	) * ddx2 &
433	+ ( &
434	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
435	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
436	) * ddy2 &
437	+ ( &
438	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
439	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
440	) * ddzw(k) &
441	- dissipation
442
443	!
444	!-- Store dissipation if needed for calculating the sgs
445	!-- particle velocities
446	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
447	turbulence ) THEN
448	diss(k,j,i) = dissipation
449	ENDIF
450
451	ENDIF
452
453	ENDDO
454	ENDDO
455	ENDDO
456	!$acc end kernels
457
458	ENDIF
459
460	!
461	!-- Boundary condition for dissipation
462	IF ( use_sgs_for_particles .OR. wang_kernel .OR. turbulence ) THEN
463	!$acc kernels present( diss, nzb_s_inner )
464	DO i = i_left, i_right
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_single_reference_value ) 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 @ 1257

Download in other formats: