Home

Context Navigation

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

Last change on this file since 790 was 790, checked in by raasch, 13 years ago
Bugfix for output of mean particle radius + preliminary works for implementing the Wang collision kernel
Property svn:keywords set to `Id`
File size: 15.3 KB

Line
1	MODULE diffusion_e_mod
2
3	!------------------------------------------------------------------------------!
4	! Current revisions:
5	! -----------------
6	! diss is also calculated in case that the Wang kernel is used
7	!
8	! Former revisions:
9	! -----------------
10	! $Id: diffusion_e.f90 790 2011-11-29 03:11:20Z raasch $
11	!
12	! 667 2010-12-23 12:06:00Z suehring/gryschka
13	! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
14	!
15	! 97 2007-06-21 08:23:15Z raasch
16	! Adjustment of mixing length calculation for the ocean version. zw added to
17	! argument list.
18	! This is also a bugfix, because the height above the topography is now
19	! used instead of the height above level k=0.
20	! theta renamed var, dpt_dz renamed dvar_dz, +new argument var_reference
21	! use_pt_reference renamed use_reference
22	!
23	! 65 2007-03-13 12:11:43Z raasch
24	! Reference temperature pt_reference can be used in buoyancy term
25	!
26	! 20 2007-02-26 00:12:32Z raasch
27	! Bugfix: ddzw dimensioned 1:nzt"+1"
28	! Calculation extended for gridpoint nzt
29	!
30	! RCS Log replace by Id keyword, revision history cleaned up
31	!
32	! Revision 1.18 2006/08/04 14:29:43 raasch
33	! dissipation is stored in extra array diss if needed later on for calculating
34	! the sgs particle velocities
35	!
36	! Revision 1.1 1997/09/19 07:40:24 raasch
37	! Initial revision
38	!
39	!
40	! Description:
41	! ------------
42	! Diffusion- and dissipation terms for the TKE
43	!------------------------------------------------------------------------------!
44
45	PRIVATE
46	PUBLIC diffusion_e
47
48
49	INTERFACE diffusion_e
50	MODULE PROCEDURE diffusion_e
51	MODULE PROCEDURE diffusion_e_ij
52	END INTERFACE diffusion_e
53
54	CONTAINS
55
56
57	!------------------------------------------------------------------------------!
58	! Call for all grid points
59	!------------------------------------------------------------------------------!
60	SUBROUTINE diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, var, &
61	var_reference, rif, tend, zu, zw )
62
63	USE control_parameters
64	USE grid_variables
65	USE indices
66	USE particle_attributes
67
68	IMPLICIT NONE
69
70	INTEGER :: i, j, k
71	REAL :: dvar_dz, l_stable, phi_m, var_reference
72	REAL :: ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt+1), &
73	l_grid(1:nzt), zu(nzb:nzt+1), zw(nzb:nzt+1)
74	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: diss, tend
75	REAL, DIMENSION(:,:), POINTER :: rif
76	REAL, DIMENSION(:,:,:), POINTER :: e, km, var
77	REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: dissipation, l, ll
78
79
80	!
81	!-- This if clause must be outside the k-loop because otherwise
82	!-- runtime errors occur with -C hopt on NEC
83	IF ( use_reference ) THEN
84
85	DO i = nxl, nxr
86	DO j = nys, nyn
87	!
88	!-- First, calculate phi-function for eventually adjusting the &
89	!-- mixing length to the prandtl mixing length
90	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
91	IF ( rif(j,i) >= 0.0 ) THEN
92	phi_m = 1.0 + 5.0 * rif(j,i)
93	ELSE
94	phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
95	ENDIF
96	ENDIF
97
98	DO k = nzb_s_inner(j,i)+1, nzt
99	!
100	!-- Calculate the mixing length (for dissipation)
101	dvar_dz = atmos_ocean_sign * &
102	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
103	IF ( dvar_dz > 0.0 ) THEN
104	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
105	SQRT( g / var_reference * dvar_dz ) + 1E-5
106	ELSE
107	l_stable = l_grid(k)
108	ENDIF
109	!
110	!-- Adjustment of the mixing length
111	IF ( wall_adjustment ) THEN
112	l(k,j) = MIN( wall_adjustment_factor * &
113	( zu(k) - zw(nzb_s_inner(j,i)) ), &
114	l_grid(k), l_stable )
115	ll(k,j) = MIN( wall_adjustment_factor * &
116	( zu(k) - zw(nzb_s_inner(j,i)) ), &
117	l_grid(k) )
118	ELSE
119	l(k,j) = MIN( l_grid(k), l_stable )
120	ll(k,j) = l_grid(k)
121	ENDIF
122	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
123	l(k,j) = MIN( l(k,j), kappa * &
124	( zu(k) - zw(nzb_s_inner(j,i)) ) &
125	/ phi_m )
126	ll(k,j) = MIN( ll(k,j), kappa * &
127	( zu(k) - zw(nzb_s_inner(j,i)) ) &
128	/ phi_m )
129	ENDIF
130
131	ENDDO
132	ENDDO
133
134	!
135	!-- Calculate the tendency terms
136	DO j = nys, nyn
137	DO k = nzb_s_inner(j,i)+1, nzt
138
139	dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
140	e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
141
142	tend(k,j,i) = tend(k,j,i) &
143	+ ( &
144	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
145	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
146	) * ddx2 &
147	+ ( &
148	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
149	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
150	) * ddy2 &
151	+ ( &
152	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
153	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
154	) * ddzw(k) &
155	- dissipation(k,j)
156
157	ENDDO
158	ENDDO
159
160	!
161	!-- Store dissipation if needed for calculating the sgs particle
162	!-- velocities
163	IF ( use_sgs_for_particles .OR. wang_collision_kernel ) THEN
164	DO j = nys, nyn
165	DO k = nzb_s_inner(j,i)+1, nzt
166	diss(k,j,i) = dissipation(k,j)
167	ENDDO
168	ENDDO
169	ENDIF
170
171	ENDDO
172
173	ELSE
174
175	DO i = nxl, nxr
176	DO j = nys, nyn
177	!
178	!-- First, calculate phi-function for eventually adjusting the &
179	!-- mixing length to the prandtl mixing length
180	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
181	IF ( rif(j,i) >= 0.0 ) THEN
182	phi_m = 1.0 + 5.0 * rif(j,i)
183	ELSE
184	phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
185	ENDIF
186	ENDIF
187
188	DO k = nzb_s_inner(j,i)+1, nzt
189	!
190	!-- Calculate the mixing length (for dissipation)
191	dvar_dz = atmos_ocean_sign * &
192	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
193	IF ( dvar_dz > 0.0 ) THEN
194	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
195	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
196	ELSE
197	l_stable = l_grid(k)
198	ENDIF
199	!
200	!-- Adjustment of the mixing length
201	IF ( wall_adjustment ) THEN
202	l(k,j) = MIN( wall_adjustment_factor * &
203	( zu(k) - zw(nzb_s_inner(j,i)) ), &
204	l_grid(k), l_stable )
205	ll(k,j) = MIN( wall_adjustment_factor * &
206	( zu(k) - zw(nzb_s_inner(j,i)) ), &
207	l_grid(k) )
208	ELSE
209	l(k,j) = MIN( l_grid(k), l_stable )
210	ll(k,j) = l_grid(k)
211	ENDIF
212	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
213	l(k,j) = MIN( l(k,j), kappa * &
214	( zu(k) - zw(nzb_s_inner(j,i)) ) &
215	/ phi_m )
216	ll(k,j) = MIN( ll(k,j), kappa * &
217	( zu(k) - zw(nzb_s_inner(j,i)) ) &
218	/ phi_m )
219	ENDIF
220
221	ENDDO
222	ENDDO
223
224	!
225	!-- Calculate the tendency terms
226	DO j = nys, nyn
227	DO k = nzb_s_inner(j,i)+1, nzt
228
229	dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
230	e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
231
232	tend(k,j,i) = tend(k,j,i) &
233	+ ( &
234	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
235	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
236	) * ddx2 &
237	+ ( &
238	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
239	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
240	) * ddy2 &
241	+ ( &
242	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
243	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
244	) * ddzw(k) &
245	- dissipation(k,j)
246
247	ENDDO
248	ENDDO
249
250	!
251	!-- Store dissipation if needed for calculating the sgs particle
252	!-- velocities
253	IF ( use_sgs_for_particles .OR. wang_collision_kernel ) THEN
254	DO j = nys, nyn
255	DO k = nzb_s_inner(j,i)+1, nzt
256	diss(k,j,i) = dissipation(k,j)
257	ENDDO
258	ENDDO
259	ENDIF
260
261	ENDDO
262
263	ENDIF
264
265	!
266	!-- Boundary condition for dissipation
267	IF ( use_sgs_for_particles .OR. wang_collision_kernel ) THEN
268	DO i = nxl, nxr
269	DO j = nys, nyn
270	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
271	ENDDO
272	ENDDO
273	ENDIF
274
275	END SUBROUTINE diffusion_e
276
277
278	!------------------------------------------------------------------------------!
279	! Call for grid point i,j
280	!------------------------------------------------------------------------------!
281	SUBROUTINE diffusion_e_ij( i, j, ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
282	var, var_reference, rif, tend, zu, zw )
283
284	USE control_parameters
285	USE grid_variables
286	USE indices
287	USE particle_attributes
288
289	IMPLICIT NONE
290
291	INTEGER :: i, j, k
292	REAL :: dvar_dz, l_stable, phi_m, var_reference
293	REAL :: ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt+1), &
294	l_grid(1:nzt), zu(nzb:nzt+1), zw(nzb:nzt+1)
295	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: diss, tend
296	REAL, DIMENSION(:,:), POINTER :: rif
297	REAL, DIMENSION(:,:,:), POINTER :: e, km, var
298	REAL, DIMENSION(nzb+1:nzt) :: dissipation, l, ll
299
300
301	!
302	!-- First, calculate phi-function for eventually adjusting the mixing length
303	!-- to the prandtl mixing length
304	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
305	IF ( rif(j,i) >= 0.0 ) THEN
306	phi_m = 1.0 + 5.0 * rif(j,i)
307	ELSE
308	phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
309	ENDIF
310	ENDIF
311
312	!
313	!-- Calculate the mixing length (for dissipation)
314	DO k = nzb_s_inner(j,i)+1, nzt
315	dvar_dz = atmos_ocean_sign * &
316	( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
317	IF ( dvar_dz > 0.0 ) THEN
318	IF ( use_reference ) THEN
319	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
320	SQRT( g / var_reference * dvar_dz ) + 1E-5
321	ELSE
322	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
323	SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5
324	ENDIF
325	ELSE
326	l_stable = l_grid(k)
327	ENDIF
328	!
329	!-- Adjustment of the mixing length
330	IF ( wall_adjustment ) THEN
331	l(k) = MIN( wall_adjustment_factor * &
332	( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k), &
333	l_stable )
334	ll(k) = MIN( wall_adjustment_factor * &
335	( zu(k) - zw(nzb_s_inner(j,i)) ), l_grid(k) )
336	ELSE
337	l(k) = MIN( l_grid(k), l_stable )
338	ll(k) = l_grid(k)
339	ENDIF
340	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
341	l(k) = MIN( l(k), kappa * &
342	( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )
343	ll(k) = MIN( ll(k), kappa * &
344	( zu(k) - zw(nzb_s_inner(j,i)) ) / phi_m )
345	ENDIF
346
347	!
348	!-- Calculate the tendency term
349	dissipation(k) = ( 0.19 + 0.74 * l(k) / ll(k) ) * e(k,j,i) * &
350	SQRT( e(k,j,i) ) / l(k)
351
352	tend(k,j,i) = tend(k,j,i) &
353	+ ( &
354	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
355	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
356	) * ddx2 &
357	+ ( &
358	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
359	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
360	) * ddy2 &
361	+ ( &
362	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
363	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
364	) * ddzw(k) &
365	- dissipation(k)
366
367	ENDDO
368
369	!
370	!-- Store dissipation if needed for calculating the sgs particle velocities
371	IF ( use_sgs_for_particles .OR. wang_collision_kernel ) THEN
372	DO k = nzb_s_inner(j,i)+1, nzt
373	diss(k,j,i) = dissipation(k)
374	ENDDO
375	!
376	!-- Boundary condition for dissipation
377	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
378	ENDIF
379
380	END SUBROUTINE diffusion_e_ij
381
382	END MODULE diffusion_e_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |