Home

Context Navigation

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

Last change on this file since 2118 was 2118, checked in by raasch, 7 years ago
all OpenACC directives and related parts removed from the code
Property svn:keywords set to `Id`
File size: 17.3 KB

Line
1	!> @file diffusion_e.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of PALM.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! 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-2017 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	! OpenACC version of subroutine removed
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: diffusion_e.f90 2118 2017-01-17 16:38:49Z raasch $
27	!
28	! 2037 2016-10-26 11:15:40Z knoop
29	! Anelastic approximation implemented
30	!
31	! 2000 2016-08-20 18:09:15Z knoop
32	! Forced header and separation lines into 80 columns
33	!
34	! 1873 2016-04-18 14:50:06Z maronga
35	! Module renamed (removed _mod)
36	!
37	! 1850 2016-04-08 13:29:27Z maronga
38	! Module renamed
39	! Adapted for modularization of microphysics
40	!
41	! 1831 2016-04-07 13:15:51Z hoffmann
42	! turbulence renamed collision_turbulence
43	!
44	! 1682 2015-10-07 23:56:08Z knoop
45	! Code annotations made doxygen readable
46	!
47	! 1374 2014-04-25 12:55:07Z raasch
48	! missing variables added to ONLY list
49	! rif removed from acc-present-list
50	!
51	! 1340 2014-03-25 19:45:13Z kanani
52	! REAL constants defined as wp-kind
53	!
54	! 1320 2014-03-20 08:40:49Z raasch
55	! ONLY-attribute added to USE-statements,
56	! kind-parameters added to all INTEGER and REAL declaration statements,
57	! kinds are defined in new module kinds,
58	! revision history before 2012 removed,
59	! comment fields (!:) to be used for variable explanations added to
60	! all variable declaration statements
61	!
62	! 1257 2013-11-08 15:18:40Z raasch
63	! openacc loop and loop vector clauses removed
64	!
65	! 1179 2013-06-14 05:57:58Z raasch
66	! use_reference renamed use_single_reference_value
67	!
68	! 1171 2013-05-30 11:27:45Z raasch
69	! use_reference-case activated in accelerator version
70	!
71	! 1128 2013-04-12 06:19:32Z raasch
72	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
73	! j_north
74	!
75	! 1065 2012-11-22 17:42:36Z hoffmann
76	! Enabled the claculation of diss in case of turbulence = .TRUE. (parameterized
77	! effects of turbulence on autoconversion and accretion in two-moments cloud
78	! physics scheme).
79	!
80	! 1036 2012-10-22 13:43:42Z raasch
81	! code put under GPL (PALM 3.9)
82	!
83	! 1015 2012-09-27 09:23:24Z raasch
84	! accelerator version (*_acc) added,
85	! adjustment of mixing length to the Prandtl mixing length at first grid point
86	! above ground removed
87	!
88	! 1010 2012-09-20 07:59:54Z raasch
89	! cpp switch __nopointer added for pointer free version
90	!
91	! 1001 2012-09-13 14:08:46Z raasch
92	! most arrays comunicated by module instead of parameter list
93	!
94	! 825 2012-02-19 03:03:44Z raasch
95	! wang_collision_kernel renamed wang_kernel
96	!
97	! Revision 1.1 1997/09/19 07:40:24 raasch
98	! Initial revision
99	!
100	!
101	! Description:
102	! ------------
103	!> Diffusion- and dissipation terms for the TKE
104	!------------------------------------------------------------------------------!
105	MODULE diffusion_e_mod
106
107
108	PRIVATE
109	PUBLIC diffusion_e
110
111
112	INTERFACE diffusion_e
113	MODULE PROCEDURE diffusion_e
114	MODULE PROCEDURE diffusion_e_ij
115	END INTERFACE diffusion_e
116
117	CONTAINS
118
119
120	!------------------------------------------------------------------------------!
121	! Description:
122	! ------------
123	!> Call for all grid points
124	!------------------------------------------------------------------------------!
125	SUBROUTINE diffusion_e( var, var_reference )
126
127	USE arrays_3d, &
128	ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw, &
129	drho_air, rho_air_zw
130
131	USE control_parameters, &
132	ONLY: atmos_ocean_sign, g, use_single_reference_value, &
133	wall_adjustment, wall_adjustment_factor
134
135	USE grid_variables, &
136	ONLY: ddx2, ddy2
137
138	USE indices, &
139	ONLY: nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzb_s_inner,&
140	nzt
141
142	USE kinds
143
144	USE microphysics_mod, &
145	ONLY: collision_turbulence
146
147	USE particle_attributes, &
148	ONLY: use_sgs_for_particles, wang_kernel
149
150	IMPLICIT NONE
151
152	INTEGER(iwp) :: i !<
153	INTEGER(iwp) :: j !<
154	INTEGER(iwp) :: k !<
155	REAL(wp) :: dvar_dz !<
156	REAL(wp) :: l_stable !<
157	REAL(wp) :: var_reference !<
158
159	#if defined( __nopointer )
160	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !<
161	#else
162	REAL(wp), DIMENSION(:,:,:), POINTER :: var !<
163	#endif
164	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dissipation !<
165	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: l !<
166	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: ll !<
167
168
169	!
170	!-- This if clause must be outside the k-loop because otherwise
171	!-- runtime errors occur with -C hopt on NEC
172	IF ( use_single_reference_value ) THEN
173
174	DO i = nxl, nxr
175	DO j = nys, nyn
176	DO k = nzb_s_inner(j,i)+1, nzt
177	!
178	!-- Calculate the mixing length (for dissipation)
179	dvar_dz = atmos_ocean_sign * &
180	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
181	IF ( dvar_dz > 0.0_wp ) THEN
182	l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
183	SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
184	ELSE
185	l_stable = l_grid(k)
186	ENDIF
187	!
188	!-- Adjustment of the mixing length
189	IF ( wall_adjustment ) THEN
190	l(k,j) = MIN( wall_adjustment_factor * &
191	( zu(k) - zw(nzb_s_inner(j,i)) ), &
192	l_grid(k), l_stable )
193	ll(k,j) = MIN( wall_adjustment_factor * &
194	( zu(k) - zw(nzb_s_inner(j,i)) ), &
195	l_grid(k) )
196	ELSE
197	l(k,j) = MIN( l_grid(k), l_stable )
198	ll(k,j) = l_grid(k)
199	ENDIF
200
201	ENDDO
202	ENDDO
203
204	!
205	!-- Calculate the tendency terms
206	DO j = nys, nyn
207	DO k = nzb_s_inner(j,i)+1, nzt
208
209	dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) * &
210	e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
211
212	tend(k,j,i) = tend(k,j,i) &
213	+ ( &
214	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
215	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
216	) * ddx2 &
217	+ ( &
218	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
219	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
220	) * ddy2 &
221	+ ( &
222	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
223	* rho_air_zw(k) &
224	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
225	* rho_air_zw(k-1) &
226	) * ddzw(k) * drho_air(k) &
227	- dissipation(k,j)
228
229	ENDDO
230	ENDDO
231
232	!
233	!-- Store dissipation if needed for calculating the sgs particle
234	!-- velocities
235	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
236	collision_turbulence ) THEN
237	DO j = nys, nyn
238	DO k = nzb_s_inner(j,i)+1, nzt
239	diss(k,j,i) = dissipation(k,j)
240	ENDDO
241	ENDDO
242	ENDIF
243
244	ENDDO
245
246	ELSE
247
248	DO i = nxl, nxr
249	DO j = nys, nyn
250	DO k = nzb_s_inner(j,i)+1, nzt
251	!
252	!-- Calculate the mixing length (for dissipation)
253	dvar_dz = atmos_ocean_sign * &
254	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
255	IF ( dvar_dz > 0.0_wp ) THEN
256	l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
257	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
258	ELSE
259	l_stable = l_grid(k)
260	ENDIF
261	!
262	!-- Adjustment of the mixing length
263	IF ( wall_adjustment ) THEN
264	l(k,j) = MIN( wall_adjustment_factor * &
265	( zu(k) - zw(nzb_s_inner(j,i)) ), &
266	l_grid(k), l_stable )
267	ll(k,j) = MIN( wall_adjustment_factor * &
268	( zu(k) - zw(nzb_s_inner(j,i)) ), &
269	l_grid(k) )
270	ELSE
271	l(k,j) = MIN( l_grid(k), l_stable )
272	ll(k,j) = l_grid(k)
273	ENDIF
274
275	ENDDO
276	ENDDO
277
278	!
279	!-- Calculate the tendency terms
280	DO j = nys, nyn
281	DO k = nzb_s_inner(j,i)+1, nzt
282
283	dissipation(k,j) = ( 0.19_wp + 0.74_wp * l(k,j) / ll(k,j) ) * &
284	e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
285
286	tend(k,j,i) = tend(k,j,i) &
287	+ ( &
288	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
289	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
290	) * ddx2 &
291	+ ( &
292	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
293	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
294	) * ddy2 &
295	+ ( &
296	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
297	* rho_air_zw(k) &
298	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
299	* rho_air_zw(k-1) &
300	) * ddzw(k) * drho_air(k) &
301	- dissipation(k,j)
302
303	ENDDO
304	ENDDO
305
306	!
307	!-- Store dissipation if needed for calculating the sgs particle
308	!-- velocities
309	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
310	collision_turbulence ) THEN
311	DO j = nys, nyn
312	DO k = nzb_s_inner(j,i)+1, nzt
313	diss(k,j,i) = dissipation(k,j)
314	ENDDO
315	ENDDO
316	ENDIF
317
318	ENDDO
319
320	ENDIF
321
322	!
323	!-- Boundary condition for dissipation
324	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
325	collision_turbulence ) THEN
326	DO i = nxl, nxr
327	DO j = nys, nyn
328	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
329	ENDDO
330	ENDDO
331	ENDIF
332
333	END SUBROUTINE diffusion_e
334
335
336	!------------------------------------------------------------------------------!
337	! Description:
338	! ------------
339	!> Call for grid point i,j
340	!------------------------------------------------------------------------------!
341	SUBROUTINE diffusion_e_ij( i, j, var, var_reference )
342
343	USE arrays_3d, &
344	ONLY: dd2zu, ddzu, ddzw, diss, e, km, l_grid, tend, zu, zw, &
345	drho_air, rho_air_zw
346
347	USE control_parameters, &
348	ONLY: atmos_ocean_sign, g, use_single_reference_value, &
349	wall_adjustment, wall_adjustment_factor
350
351	USE grid_variables, &
352	ONLY: ddx2, ddy2
353
354	USE indices, &
355	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzb_s_inner, nzt
356
357	USE kinds
358
359	USE microphysics_mod, &
360	ONLY: collision_turbulence
361
362	USE particle_attributes, &
363	ONLY: use_sgs_for_particles, wang_kernel
364
365	IMPLICIT NONE
366
367	INTEGER(iwp) :: i !<
368	INTEGER(iwp) :: j !<
369	INTEGER(iwp) :: k !<
370	REAL(wp) :: dvar_dz !<
371	REAL(wp) :: l_stable !<
372	REAL(wp) :: var_reference !<
373
374	#if defined( __nopointer )
375	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !<
376	#else
377	REAL(wp), DIMENSION(:,:,:), POINTER :: var !<
378	#endif
379	REAL(wp), DIMENSION(nzb+1:nzt) :: dissipation !<
380	REAL(wp), DIMENSION(nzb+1:nzt) :: l !<
381	REAL(wp), DIMENSION(nzb+1:nzt) :: ll !<
382
383
384	!
385	!-- Calculate the mixing length (for dissipation)
386	DO k = nzb_s_inner(j,i)+1, nzt
387	dvar_dz = atmos_ocean_sign * &
388	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
389	IF ( dvar_dz > 0.0_wp ) THEN
390	IF ( use_single_reference_value ) THEN
391	l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
392	SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
393	ELSE
394	l_stable = 0.76_wp * SQRT( e(k,j,i) ) / &
395	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
396	ENDIF
397	ELSE
398	l_stable = l_grid(k)
399	ENDIF
400	!
401	!-- Adjustment of the mixing length
402	IF ( wall_adjustment ) THEN
403	l(k) = MIN( wall_adjustment_factor * &
404	( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k), &
405	l_stable )
406	ll(k) = MIN( wall_adjustment_factor * &
407	( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k) )
408	ELSE
409	l(k) = MIN( l_grid(k), l_stable )
410	ll(k) = l_grid(k)
411	ENDIF
412	!
413	!-- Calculate the tendency term
414	dissipation(k) = ( 0.19_wp + 0.74_wp * l(k) / ll(k) ) * e(k,j,i) * &
415	SQRT( e(k,j,i) ) / l(k)
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	* rho_air_zw(k) &
429	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
430	* rho_air_zw(k-1) &
431	) * ddzw(k) * drho_air(k) &
432	- dissipation(k)
433
434	ENDDO
435
436	!
437	!-- Store dissipation if needed for calculating the sgs particle velocities
438	IF ( use_sgs_for_particles .OR. wang_kernel .OR. &
439	collision_turbulence ) THEN
440	DO k = nzb_s_inner(j,i)+1, nzt
441	diss(k,j,i) = dissipation(k)
442	ENDDO
443	!
444	!-- Boundary condition for dissipation
445	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
446	ENDIF
447
448	END SUBROUTINE diffusion_e_ij
449
450	END MODULE diffusion_e_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |