Home

Context Navigation

source: palm/trunk/SOURCE/diffusion_s.f90 @ 1117

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

Line
1	MODULE diffusion_s_mod
2
3	!--------------------------------------------------------------------------------!
4	! This file is part of PALM.
5	!
6	! PALM is free software: you can redistribute it and/or modify it under the terms
7	! of the GNU General Public License as published by the Free Software Foundation,
8	! either version 3 of the License, or (at your option) any later version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2012 Leibniz University Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: diffusion_s.f90 1093 2013-02-02 12:58:49Z suehring $
27	!
28	! 1092 2013-02-02 11:24:22Z raasch
29	! unused variables removed
30	!
31	! 1036 2012-10-22 13:43:42Z raasch
32	! code put under GPL (PALM 3.9)
33	!
34	! 1015 2012-09-27 09:23:24Z raasch
35	! accelerator version (*_acc) added
36	!
37	! 1010 2012-09-20 07:59:54Z raasch
38	! cpp switch __nopointer added for pointer free version
39	!
40	! 1001 2012-09-13 14:08:46Z raasch
41	! some arrays comunicated by module instead of parameter list
42	!
43	! 667 2010-12-23 12:06:00Z suehring/gryschka
44	! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
45	!
46	! 183 2008-08-04 15:39:12Z letzel
47	! bugfix: calculation of fluxes at vertical surfaces
48	!
49	! 129 2007-10-30 12:12:24Z letzel
50	! replace wall_heatflux by wall_s_flux that is now included in the parameter
51	! list, bugfix for assignment of fluxes at walls
52	!
53	! 20 2007-02-26 00:12:32Z raasch
54	! Bugfix: ddzw dimensioned 1:nzt"+1"
55	! Calculation extended for gridpoint nzt, fluxes can be given at top,
56	! +s_flux_t in parameter list, s_flux renamed s_flux_b
57	!
58	! RCS Log replace by Id keyword, revision history cleaned up
59	!
60	! Revision 1.8 2006/02/23 10:34:17 raasch
61	! nzb_2d replaced by nzb_s_outer in horizontal diffusion and by nzb_s_inner
62	! or nzb_diff_s_inner, respectively, in vertical diffusion, prescribed surface
63	! fluxes at vertically oriented topography
64	!
65	! Revision 1.1 2000/04/13 14:54:02 schroeter
66	! Initial revision
67	!
68	!
69	! Description:
70	! ------------
71	! Diffusion term of scalar quantities (temperature and water content)
72	!------------------------------------------------------------------------------!
73
74	PRIVATE
75	PUBLIC diffusion_s, diffusion_s_acc
76
77	INTERFACE diffusion_s
78	MODULE PROCEDURE diffusion_s
79	MODULE PROCEDURE diffusion_s_ij
80	END INTERFACE diffusion_s
81
82	INTERFACE diffusion_s_acc
83	MODULE PROCEDURE diffusion_s_acc
84	END INTERFACE diffusion_s_acc
85
86	CONTAINS
87
88
89	!------------------------------------------------------------------------------!
90	! Call for all grid points
91	!------------------------------------------------------------------------------!
92	SUBROUTINE diffusion_s( s, s_flux_b, s_flux_t, wall_s_flux )
93
94	USE arrays_3d
95	USE control_parameters
96	USE grid_variables
97	USE indices
98
99	IMPLICIT NONE
100
101	INTEGER :: i, j, k
102	REAL :: wall_s_flux(0:4)
103	REAL, DIMENSION(nysg:nyng,nxlg:nxrg) :: s_flux_b, s_flux_t
104	#if defined( __nopointer )
105	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: s
106	#else
107	REAL, DIMENSION(:,:,:), POINTER :: s
108	#endif
109
110	DO i = nxl, nxr
111	DO j = nys,nyn
112	!
113	!-- Compute horizontal diffusion
114	DO k = nzb_s_outer(j,i)+1, nzt
115
116	tend(k,j,i) = tend(k,j,i) &
117	+ 0.5 * ( &
118	( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
119	- ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
120	) * ddx2 &
121	+ 0.5 * ( &
122	( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
123	- ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
124	) * ddy2
125	ENDDO
126
127	!
128	!-- Apply prescribed horizontal wall heatflux where necessary
129	IF ( ( wall_w_x(j,i) .NE. 0.0 ) .OR. ( wall_w_y(j,i) .NE. 0.0 ) ) &
130	THEN
131	DO k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i)
132
133	tend(k,j,i) = tend(k,j,i) &
134	+ ( fwxp(j,i) * 0.5 * &
135	( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
136	+ ( 1.0 - fwxp(j,i) ) * wall_s_flux(1) &
137	-fwxm(j,i) * 0.5 * &
138	( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
139	+ ( 1.0 - fwxm(j,i) ) * wall_s_flux(2) &
140	) * ddx2 &
141	+ ( fwyp(j,i) * 0.5 * &
142	( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
143	+ ( 1.0 - fwyp(j,i) ) * wall_s_flux(3) &
144	-fwym(j,i) * 0.5 * &
145	( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
146	+ ( 1.0 - fwym(j,i) ) * wall_s_flux(4) &
147	) * ddy2
148	ENDDO
149	ENDIF
150
151	!
152	!-- Compute vertical diffusion. In case that surface fluxes have been
153	!-- prescribed or computed at bottom and/or top, index k starts/ends at
154	!-- nzb+2 or nzt-1, respectively.
155	DO k = nzb_diff_s_inner(j,i), nzt_diff
156
157	tend(k,j,i) = tend(k,j,i) &
158	+ 0.5 * ( &
159	( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
160	- ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k) &
161	) * ddzw(k)
162	ENDDO
163
164	!
165	!-- Vertical diffusion at the first computational gridpoint along
166	!-- z-direction
167	IF ( use_surface_fluxes ) THEN
168
169	k = nzb_s_inner(j,i)+1
170
171	tend(k,j,i) = tend(k,j,i) &
172	+ ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) ) &
173	* ( s(k+1,j,i)-s(k,j,i) ) &
174	* ddzu(k+1) &
175	+ s_flux_b(j,i) &
176	) * ddzw(k)
177
178	ENDIF
179
180	!
181	!-- Vertical diffusion at the last computational gridpoint along
182	!-- z-direction
183	IF ( use_top_fluxes ) THEN
184
185	k = nzt
186
187	tend(k,j,i) = tend(k,j,i) &
188	+ ( - s_flux_t(j,i) &
189	- 0.5 * ( kh(k-1,j,i)+kh(k,j,i) ) &
190	* ( s(k,j,i)-s(k-1,j,i) ) &
191	* ddzu(k) &
192	) * ddzw(k)
193
194	ENDIF
195
196	ENDDO
197	ENDDO
198
199	END SUBROUTINE diffusion_s
200
201
202	!------------------------------------------------------------------------------!
203	! Call for all grid points - accelerator version
204	!------------------------------------------------------------------------------!
205	SUBROUTINE diffusion_s_acc( s, s_flux_b, s_flux_t, wall_s_flux )
206
207	USE arrays_3d
208	USE control_parameters
209	USE grid_variables
210	USE indices
211
212	IMPLICIT NONE
213
214	INTEGER :: i, j, k
215	REAL :: wall_s_flux(0:4)
216	REAL, DIMENSION(nysg:nyng,nxlg:nxrg) :: s_flux_b, s_flux_t
217	#if defined( __nopointer )
218	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: s
219	#else
220	REAL, DIMENSION(:,:,:), POINTER :: s
221	#endif
222
223	!$acc kernels present( ddzu, ddzw, fwxm, fwxp, fwym, fwyp, kh ) &
224	!$acc present( nzb_diff_s_inner, nzb_s_inner, nzb_s_outer, s ) &
225	!$acc present( s_flux_b, s_flux_t, tend, wall_s_flux ) &
226	!$acc present( wall_w_x, wall_w_y )
227	!$acc loop
228	DO i = nxl, nxr
229	DO j = nys,nyn
230	!
231	!-- Compute horizontal diffusion
232	!$acc loop vector( 32 )
233	DO k = 1, nzt
234	IF ( k > nzb_s_outer(j,i) ) THEN
235
236	tend(k,j,i) = tend(k,j,i) &
237	+ 0.5 * ( &
238	( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
239	- ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
240	) * ddx2 &
241	+ 0.5 * ( &
242	( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
243	- ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
244	) * ddy2
245	ENDIF
246	ENDDO
247
248	!
249	!-- Apply prescribed horizontal wall heatflux where necessary
250	!$acc loop vector(32)
251	DO k = 1, nzt
252	IF ( k > nzb_s_inner(j,i) .AND. k <= nzb_s_outer(j,i) .AND. &
253	( wall_w_x(j,i) /= 0.0 .OR. wall_w_y(j,i) /= 0.0 ) ) &
254	THEN
255	tend(k,j,i) = tend(k,j,i) &
256	+ ( fwxp(j,i) * 0.5 * &
257	( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
258	+ ( 1.0 - fwxp(j,i) ) * wall_s_flux(1) &
259	-fwxm(j,i) * 0.5 * &
260	( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
261	+ ( 1.0 - fwxm(j,i) ) * wall_s_flux(2) &
262	) * ddx2 &
263	+ ( fwyp(j,i) * 0.5 * &
264	( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
265	+ ( 1.0 - fwyp(j,i) ) * wall_s_flux(3) &
266	-fwym(j,i) * 0.5 * &
267	( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
268	+ ( 1.0 - fwym(j,i) ) * wall_s_flux(4) &
269	) * ddy2
270	ENDIF
271	ENDDO
272
273	!
274	!-- Compute vertical diffusion. In case that surface fluxes have been
275	!-- prescribed or computed at bottom and/or top, index k starts/ends at
276	!-- nzb+2 or nzt-1, respectively.
277	!$acc loop vector( 32 )
278	DO k = 1, nzt_diff
279	IF ( k >= nzb_diff_s_inner(j,i) ) THEN
280	tend(k,j,i) = tend(k,j,i) &
281	+ 0.5 * ( &
282	( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
283	- ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k) &
284	) * ddzw(k)
285	ENDIF
286	ENDDO
287
288	!
289	!-- Vertical diffusion at the first computational gridpoint along
290	!-- z-direction
291	!$acc loop vector( 32 )
292	DO k = 1, nzt
293	IF ( use_surface_fluxes .AND. k == nzb_s_inner(j,i)+1 ) THEN
294	tend(k,j,i) = tend(k,j,i) &
295	+ ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) ) &
296	* ( s(k+1,j,i)-s(k,j,i) ) &
297	* ddzu(k+1) &
298	+ s_flux_b(j,i) &
299	) * ddzw(k)
300	ENDIF
301
302	!
303	!-- Vertical diffusion at the last computational gridpoint along
304	!-- z-direction
305	IF ( use_top_fluxes .AND. k == nzt ) THEN
306	tend(k,j,i) = tend(k,j,i) &
307	+ ( - s_flux_t(j,i) &
308	- 0.5 * ( kh(k-1,j,i)+kh(k,j,i) )&
309	* ( s(k,j,i)-s(k-1,j,i) ) &
310	* ddzu(k) &
311	) * ddzw(k)
312	ENDIF
313	ENDDO
314
315	ENDDO
316	ENDDO
317	!$acc end kernels
318
319	END SUBROUTINE diffusion_s_acc
320
321
322	!------------------------------------------------------------------------------!
323	! Call for grid point i,j
324	!------------------------------------------------------------------------------!
325	SUBROUTINE diffusion_s_ij( i, j, s, s_flux_b, s_flux_t, wall_s_flux )
326
327	USE arrays_3d
328	USE control_parameters
329	USE grid_variables
330	USE indices
331
332	IMPLICIT NONE
333
334	INTEGER :: i, j, k
335	REAL :: wall_s_flux(0:4)
336	REAL, DIMENSION(nysg:nyng,nxlg:nxrg) :: s_flux_b, s_flux_t
337	#if defined( __nopointer )
338	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: s
339	#else
340	REAL, DIMENSION(:,:,:), POINTER :: s
341	#endif
342
343	!
344	!-- Compute horizontal diffusion
345	DO k = nzb_s_outer(j,i)+1, nzt
346
347	tend(k,j,i) = tend(k,j,i) &
348	+ 0.5 * ( &
349	( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
350	- ( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
351	) * ddx2 &
352	+ 0.5 * ( &
353	( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
354	- ( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
355	) * ddy2
356	ENDDO
357
358	!
359	!-- Apply prescribed horizontal wall heatflux where necessary
360	IF ( ( wall_w_x(j,i) .NE. 0.0 ) .OR. ( wall_w_y(j,i) .NE. 0.0 ) ) &
361	THEN
362	DO k = nzb_s_inner(j,i)+1, nzb_s_outer(j,i)
363
364	tend(k,j,i) = tend(k,j,i) &
365	+ ( fwxp(j,i) * 0.5 * &
366	( kh(k,j,i) + kh(k,j,i+1) ) * ( s(k,j,i+1)-s(k,j,i) ) &
367	+ ( 1.0 - fwxp(j,i) ) * wall_s_flux(1) &
368	-fwxm(j,i) * 0.5 * &
369	( kh(k,j,i) + kh(k,j,i-1) ) * ( s(k,j,i)-s(k,j,i-1) ) &
370	+ ( 1.0 - fwxm(j,i) ) * wall_s_flux(2) &
371	) * ddx2 &
372	+ ( fwyp(j,i) * 0.5 * &
373	( kh(k,j,i) + kh(k,j+1,i) ) * ( s(k,j+1,i)-s(k,j,i) ) &
374	+ ( 1.0 - fwyp(j,i) ) * wall_s_flux(3) &
375	-fwym(j,i) * 0.5 * &
376	( kh(k,j,i) + kh(k,j-1,i) ) * ( s(k,j,i)-s(k,j-1,i) ) &
377	+ ( 1.0 - fwym(j,i) ) * wall_s_flux(4) &
378	) * ddy2
379	ENDDO
380	ENDIF
381
382	!
383	!-- Compute vertical diffusion. In case that surface fluxes have been
384	!-- prescribed or computed at bottom and/or top, index k starts/ends at
385	!-- nzb+2 or nzt-1, respectively.
386	DO k = nzb_diff_s_inner(j,i), nzt_diff
387
388	tend(k,j,i) = tend(k,j,i) &
389	+ 0.5 * ( &
390	( kh(k,j,i) + kh(k+1,j,i) ) * ( s(k+1,j,i)-s(k,j,i) ) * ddzu(k+1) &
391	- ( kh(k,j,i) + kh(k-1,j,i) ) * ( s(k,j,i)-s(k-1,j,i) ) * ddzu(k) &
392	) * ddzw(k)
393	ENDDO
394
395	!
396	!-- Vertical diffusion at the first computational gridpoint along z-direction
397	IF ( use_surface_fluxes ) THEN
398
399	k = nzb_s_inner(j,i)+1
400
401	tend(k,j,i) = tend(k,j,i) + ( 0.5 * ( kh(k,j,i)+kh(k+1,j,i) ) &
402	* ( s(k+1,j,i)-s(k,j,i) ) &
403	* ddzu(k+1) &
404	+ s_flux_b(j,i) &
405	) * ddzw(k)
406
407	ENDIF
408
409	!
410	!-- Vertical diffusion at the last computational gridpoint along z-direction
411	IF ( use_top_fluxes ) THEN
412
413	k = nzt
414
415	tend(k,j,i) = tend(k,j,i) + ( - s_flux_t(j,i) &
416	- 0.5 * ( kh(k-1,j,i)+kh(k,j,i) ) &
417	* ( s(k,j,i)-s(k-1,j,i) ) &
418	* ddzu(k) &
419	) * ddzw(k)
420
421	ENDIF
422
423	END SUBROUTINE diffusion_s_ij
424
425	END MODULE diffusion_s_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |