Home

Context Navigation

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

Last change on this file since 17 was 4, checked in by raasch, 18 years ago
Id keyword set as property for all *.f90 files
Property svn:keywords set to `Id`
File size: 9.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |