Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |