Home

Context Navigation

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

Last change on this file since 1011 was 1011, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 15.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |