Home

Context Navigation

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

Last change on this file since 1001 was 1001, checked in by raasch, 12 years ago
leapfrog timestep scheme and upstream-spline advection scheme completely removed from the code, reading of dt_fixed from restart file removed
Property svn:keywords set to `Id`
File size: 14.8 KB

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