Home

Context Navigation

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

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