Home

Context Navigation

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

Last change on this file since 1010 was 1010, checked in by raasch, 12 years ago
pointer free version can be generated with cpp switch nopointer
Property svn:keywords set to `Id`
File size: 15.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |