Home

Context Navigation

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

Last change on this file since 57 was 57, checked in by raasch, 17 years ago
preliminary update of further changes, advec_particles is not running!
Property svn:keywords set to `Id`
File size: 9.7 KB

Line
1	MODULE diffusion_e_mod
2
3	!------------------------------------------------------------------------------!
4	! Actual revisions:
5	! -----------------
6	! Reference temperature pt_reference can be used in buoyancy term
7	!
8	! Former revisions:
9	! -----------------
10	! $Id: diffusion_e.f90 57 2007-03-09 12:05:41Z raasch $
11	!
12	! 20 2007-02-26 00:12:32Z raasch
13	! Bugfix: ddzw dimensioned 1:nzt"+1"
14	! Calculation extended for gridpoint nzt
15	!
16	! RCS Log replace by Id keyword, revision history cleaned up
17	!
18	! Revision 1.18 2006/08/04 14:29:43 raasch
19	! dissipation is stored in extra array diss if needed later on for calculating
20	! the sgs particle velocities
21	!
22	! Revision 1.1 1997/09/19 07:40:24 raasch
23	! Initial revision
24	!
25	!
26	! Description:
27	! ------------
28	! Diffusion- and dissipation terms for the TKE
29	!------------------------------------------------------------------------------!
30
31	PRIVATE
32	PUBLIC diffusion_e
33
34
35	INTERFACE diffusion_e
36	MODULE PROCEDURE diffusion_e
37	MODULE PROCEDURE diffusion_e_ij
38	END INTERFACE diffusion_e
39
40	CONTAINS
41
42
43	!------------------------------------------------------------------------------!
44	! Call for all grid points
45	!------------------------------------------------------------------------------!
46	SUBROUTINE diffusion_e( ddzu, dd2zu, ddzw, diss, e, km, l_grid, theta, &
47	rif, tend, zu )
48
49	USE control_parameters
50	USE grid_variables
51	USE indices
52	USE particle_attributes
53
54	IMPLICIT NONE
55
56	INTEGER :: i, j, k
57	REAL :: dpt_dz, l_stable, phi_m
58	REAL :: ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt+1), &
59	l_grid(1:nzt), zu(0:nzt+1)
60	REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) :: diss, tend
61	REAL, DIMENSION(:,:), POINTER :: rif
62	REAL, DIMENSION(:,:,:), POINTER :: e, km, theta
63	REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: dissipation, l, ll
64
65
66	DO i = nxl, nxr
67	DO j = nys, nyn
68	!
69	!-- First, calculate phi-function for eventually adjusting the &
70	!-- mixing length to the prandtl mixing length
71	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
72	IF ( rif(j,i) >= 0.0 ) THEN
73	phi_m = 1.0 + 5.0 * rif(j,i)
74	ELSE
75	phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
76	ENDIF
77	ENDIF
78
79	DO k = nzb_s_inner(j,i)+1, nzt
80	!
81	!-- Calculate the mixing length (for dissipation)
82	dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
83	IF ( dpt_dz > 0.0 ) THEN
84	IF ( use_pt_reference ) THEN
85	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
86	SQRT( g / pt_reference * dpt_dz ) + 1E-5
87	ELSE
88	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
89	SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
90	ENDIF
91	ELSE
92	l_stable = l_grid(k)
93	ENDIF
94	!
95	!-- Adjustment of the mixing length
96	IF ( wall_adjustment ) THEN
97	l(k,j) = MIN( wall_adjustment_factor * zu(k), l_grid(k), &
98	l_stable )
99	ll(k,j) = MIN( wall_adjustment_factor * zu(k), l_grid(k) )
100	ELSE
101	l(k,j) = MIN( l_grid(k), l_stable )
102	ll(k,j) = l_grid(k)
103	ENDIF
104	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
105	l(k,j) = MIN( l(k,j), kappa * zu(k) / phi_m )
106	ll(k,j) = MIN( ll(k,j), kappa * zu(k) / phi_m )
107	ENDIF
108
109	ENDDO
110	ENDDO
111	!
112	!-- Calculate the tendency terms
113	DO j = nys, nyn
114	DO k = nzb_s_inner(j,i)+1, nzt
115
116	dissipation(k,j) = ( 0.19 + 0.74 * l(k,j) / ll(k,j) ) * &
117	e(k,j,i) * SQRT( e(k,j,i) ) / l(k,j)
118
119	tend(k,j,i) = tend(k,j,i) &
120	+ ( &
121	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
122	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
123	) * ddx2 &
124	+ ( &
125	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
126	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
127	) * ddy2 &
128	+ ( &
129	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
130	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
131	) * ddzw(k) &
132	- dissipation(k,j)
133
134	ENDDO
135	ENDDO
136
137	!
138	!-- Store dissipation if needed for calculating the sgs particle
139	!-- velocities
140	IF ( use_sgs_for_particles ) THEN
141	DO j = nys, nyn
142	DO k = nzb_s_inner(j,i)+1, nzt
143	diss(k,j,i) = dissipation(k,j)
144	ENDDO
145	ENDDO
146	ENDIF
147
148	ENDDO
149
150	!
151	!-- Boundary condition for dissipation
152	IF ( use_sgs_for_particles ) THEN
153	DO i = nxl, nxr
154	DO j = nys, nyn
155	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
156	ENDDO
157	ENDDO
158	ENDIF
159
160	END SUBROUTINE diffusion_e
161
162
163	!------------------------------------------------------------------------------!
164	! Call for grid point i,j
165	!------------------------------------------------------------------------------!
166	SUBROUTINE diffusion_e_ij( i, j, ddzu, dd2zu, ddzw, diss, e, km, l_grid, &
167	theta, rif, tend, zu )
168
169	USE control_parameters
170	USE grid_variables
171	USE indices
172	USE particle_attributes
173
174	IMPLICIT NONE
175
176	INTEGER :: i, j, k
177	REAL :: dpt_dz, l_stable, phi_m
178	REAL :: ddzu(1:nzt+1), dd2zu(1:nzt), ddzw(1:nzt+1), &
179	l_grid(1:nzt), zu(0:nzt+1)
180	REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) :: diss, tend
181	REAL, DIMENSION(:,:), POINTER :: rif
182	REAL, DIMENSION(:,:,:), POINTER :: e, km, theta
183	REAL, DIMENSION(nzb+1:nzt) :: dissipation, l, ll
184
185
186	!
187	!-- First, calculate phi-function for eventually adjusting the mixing length
188	!-- to the prandtl mixing length
189	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
190	IF ( rif(j,i) >= 0.0 ) THEN
191	phi_m = 1.0 + 5.0 * rif(j,i)
192	ELSE
193	phi_m = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
194	ENDIF
195	ENDIF
196
197	!
198	!-- Calculate the mixing length (for dissipation)
199	DO k = nzb_s_inner(j,i)+1, nzt
200	dpt_dz = ( theta(k+1,j,i) - theta(k-1,j,i) ) * dd2zu(k)
201	IF ( dpt_dz > 0.0 ) THEN
202	IF ( use_pt_reference ) THEN
203	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
204	SQRT( g / pt_reference * dpt_dz ) + 1E-5
205	ELSE
206	l_stable = 0.76 * SQRT( e(k,j,i) ) / &
207	SQRT( g / theta(k,j,i) * dpt_dz ) + 1E-5
208	ENDIF
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) = MIN( wall_adjustment_factor * zu(k), l_grid(k), l_stable )
216	ll(k) = MIN( wall_adjustment_factor * zu(k), l_grid(k) )
217	ELSE
218	l(k) = MIN( l_grid(k), l_stable )
219	ll(k) = l_grid(k)
220	ENDIF
221	IF ( adjust_mixing_length .AND. prandtl_layer ) THEN
222	l(k) = MIN( l(k), kappa * zu(k) / phi_m )
223	ll(k) = MIN( ll(k), kappa * zu(k) / phi_m )
224	ENDIF
225
226	!
227	!-- Calculate the tendency term
228	dissipation(k) = ( 0.19 + 0.74 * l(k) / ll(k) ) * e(k,j,i) * &
229	SQRT( e(k,j,i) ) / l(k)
230
231	tend(k,j,i) = tend(k,j,i) &
232	+ ( &
233	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
234	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
235	) * ddx2 &
236	+ ( &
237	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
238	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
239	) * ddy2 &
240	+ ( &
241	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
242	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
243	) * ddzw(k) &
244	- dissipation(k)
245
246	ENDDO
247
248	!
249	!-- Store dissipation if needed for calculating the sgs particle velocities
250	IF ( use_sgs_for_particles ) THEN
251	DO k = nzb_s_inner(j,i)+1, nzt
252	diss(k,j,i) = dissipation(k)
253	ENDDO
254	!
255	!-- Boundary condition for dissipation
256	diss(nzb_s_inner(j,i),j,i) = diss(nzb_s_inner(j,i)+1,j,i)
257	ENDIF
258
259	END SUBROUTINE diffusion_e_ij
260
261	END MODULE diffusion_e_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |