Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |