Home

Context Navigation

source: palm/trunk/SOURCE/buoyancy.f90 @ 1433

Last change on this file since 1433 was 1375, checked in by raasch, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 13.5 KB

Line
1	MODULE buoyancy_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-2014 Leibniz Universitaet Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! ------------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: buoyancy.f90 1375 2014-04-25 13:07:08Z keck $
27	!
28	! 1374 2014-04-25 12:55:07Z raasch
29	! missing variables added to ONLY list
30	!
31	! 1365 2014-04-22 15:03:56Z boeske
32	! Calculation of reference state in subroutine calc_mean_profile moved to
33	! subroutine time_integration,
34	! subroutine calc_mean_profile moved to new file calc_mean_profile.f90
35	!
36	! 1353 2014-04-08 15:21:23Z heinze
37	! REAL constants provided with KIND-attribute
38	!
39	! 1320 2014-03-20 08:40:49Z raasch
40	! ONLY-attribute added to USE-statements,
41	! kind-parameters added to all INTEGER and REAL declaration statements,
42	! kinds are defined in new module kinds,
43	! revision history before 2012 removed,
44	! comment fields (!:) to be used for variable explanations added to
45	! all variable declaration statements
46	!
47	! 1257 2013-11-08 15:18:40Z raasch
48	! vector length (32) removed from openacc clause
49	!
50	! 1241 2013-10-30 11:36:58Z heinze
51	! Generalize calc_mean_profile for wider use: use additional steering
52	! character loc
53	!
54	! 1179 2013-06-14 05:57:58Z raasch
55	! steering of reference state revised (var_reference and pr removed from
56	! parameter list), use_reference renamed use_single_reference_value
57	!
58	! 1171 2013-05-30 11:27:45Z raasch
59	! openacc statements added to use_reference-case in accelerator version
60	!
61	! 1153 2013-05-10 14:33:08Z raasch
62	! code adjustments of accelerator version required by PGI 12.3 / CUDA 5.0
63	!
64	! 1128 2013-04-12 06:19:32Z raasch
65	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
66	! j_north
67	!
68	! 1036 2012-10-22 13:43:42Z raasch
69	! code put under GPL (PALM 3.9)
70	!
71	! 1015 2012-09-27 09:23:24Z raasch
72	! accelerator version (*_acc) added
73	!
74	! 1010 2012-09-20 07:59:54Z raasch
75	! cpp switch __nopointer added for pointer free version
76	!
77	! Revision 1.1 1997/08/29 08:56:48 raasch
78	! Initial revision
79	!
80	!
81	! Description:
82	! ------------
83	! Buoyancy term of the third component of the equation of motion.
84	! WARNING: humidity is not regarded when using a sloping surface!
85	!------------------------------------------------------------------------------!
86
87	PRIVATE
88	PUBLIC buoyancy, buoyancy_acc
89
90	INTERFACE buoyancy
91	MODULE PROCEDURE buoyancy
92	MODULE PROCEDURE buoyancy_ij
93	END INTERFACE buoyancy
94
95	INTERFACE buoyancy_acc
96	MODULE PROCEDURE buoyancy_acc
97	END INTERFACE buoyancy_acc
98
99	CONTAINS
100
101
102	!------------------------------------------------------------------------------!
103	! Call for all grid points
104	!------------------------------------------------------------------------------!
105	SUBROUTINE buoyancy( var, wind_component )
106
107	USE arrays_3d, &
108	ONLY: pt, pt_slope_ref, ref_state, tend
109
110	USE control_parameters, &
111	ONLY: atmos_ocean_sign, cos_alpha_surface, g, message_string, &
112	pt_surface, sin_alpha_surface, sloping_surface
113
114	USE indices, &
115	ONLY: nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nzb, &
116	nzb_s_inner, nzt
117
118	USE kinds
119
120	USE pegrid
121
122
123	IMPLICIT NONE
124
125	INTEGER(iwp) :: i !:
126	INTEGER(iwp) :: j !:
127	INTEGER(iwp) :: k !:
128	INTEGER(iwp) :: wind_component !:
129
130	#if defined( __nopointer )
131	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !:
132	#else
133	REAL(wp), DIMENSION(:,:,:), POINTER :: var
134	#endif
135
136
137	IF ( .NOT. sloping_surface ) THEN
138	!
139	!-- Normal case: horizontal surface
140	DO i = nxl, nxr
141	DO j = nys, nyn
142	DO k = nzb_s_inner(j,i)+1, nzt-1
143	tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5_wp * &
144	( &
145	( var(k,j,i) - ref_state(k) ) / ref_state(k) + &
146	( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1) &
147	)
148	ENDDO
149	ENDDO
150	ENDDO
151
152	ELSE
153	!
154	!-- Buoyancy term for a surface with a slope in x-direction. The equations
155	!-- for both the u and w velocity-component contain proportionate terms.
156	!-- Temperature field at time t=0 serves as environmental temperature.
157	!-- Reference temperature (pt_surface) is the one at the lower left corner
158	!-- of the total domain.
159	IF ( wind_component == 1 ) THEN
160
161	DO i = nxlu, nxr
162	DO j = nys, nyn
163	DO k = nzb_s_inner(j,i)+1, nzt-1
164	tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface * &
165	0.5_wp * ( ( pt(k,j,i-1) + pt(k,j,i) ) &
166	- ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
167	) / pt_surface
168	ENDDO
169	ENDDO
170	ENDDO
171
172	ELSEIF ( wind_component == 3 ) THEN
173
174	DO i = nxl, nxr
175	DO j = nys, nyn
176	DO k = nzb_s_inner(j,i)+1, nzt-1
177	tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface * &
178	0.5_wp * ( ( pt(k,j,i) + pt(k+1,j,i) ) &
179	- ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
180	) / pt_surface
181	ENDDO
182	ENDDO
183	ENDDO
184
185	ELSE
186
187	WRITE( message_string, * ) 'no term for component "', &
188	wind_component,'"'
189	CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
190
191	ENDIF
192
193	ENDIF
194
195	END SUBROUTINE buoyancy
196
197
198	!------------------------------------------------------------------------------!
199	! Call for all grid points - accelerator version
200	!------------------------------------------------------------------------------!
201	SUBROUTINE buoyancy_acc( var, wind_component )
202
203	USE arrays_3d, &
204	ONLY: pt, pt_slope_ref, ref_state, tend
205
206	USE control_parameters, &
207	ONLY: atmos_ocean_sign, cos_alpha_surface, g, message_string, &
208	pt_surface, sin_alpha_surface, sloping_surface
209
210	USE indices, &
211	ONLY: i_left, i_right, j_north, j_south, nxl, nxlg, nxlu, nxr, &
212	nxrg, nyn, nyng, nys, nysg, nzb, nzb_s_inner, nzt
213
214	USE kinds
215
216	USE pegrid
217
218
219	IMPLICIT NONE
220
221	INTEGER(iwp) :: i !:
222	INTEGER(iwp) :: j !:
223	INTEGER(iwp) :: k !:
224	INTEGER(iwp) :: wind_component !:
225
226	#if defined( __nopointer )
227	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !:
228	#else
229	REAL(wp), DIMENSION(:,:,:), POINTER :: var
230	#endif
231
232
233	IF ( .NOT. sloping_surface ) THEN
234	!
235	!-- Normal case: horizontal surface
236	!$acc kernels present( nzb_s_inner, ref_state, tend, var )
237	!$acc loop
238	DO i = i_left, i_right
239	DO j = j_south, j_north
240	!$acc loop independent vector
241	DO k = nzb_s_inner(j,i)+1, nzt-1
242	tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5_wp * &
243	( &
244	( var(k,j,i) - ref_state(k) ) / ref_state(k) + &
245	( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1) &
246	)
247	ENDDO
248	ENDDO
249	ENDDO
250	!$acc end kernels
251
252	ELSE
253	!
254	!-- Buoyancy term for a surface with a slope in x-direction. The equations
255	!-- for both the u and w velocity-component contain proportionate terms.
256	!-- Temperature field at time t=0 serves as environmental temperature.
257	!-- Reference temperature (pt_surface) is the one at the lower left corner
258	!-- of the total domain.
259	IF ( wind_component == 1 ) THEN
260
261	DO i = nxlu, nxr
262	DO j = nys, nyn
263	DO k = nzb_s_inner(j,i)+1, nzt-1
264	tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface * &
265	0.5_wp * ( ( pt(k,j,i-1) + pt(k,j,i) ) &
266	- ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
267	) / pt_surface
268	ENDDO
269	ENDDO
270	ENDDO
271
272	ELSEIF ( wind_component == 3 ) THEN
273
274	DO i = nxl, nxr
275	DO j = nys, nyn
276	DO k = nzb_s_inner(j,i)+1, nzt-1
277	tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface * &
278	0.5_wp * ( ( pt(k,j,i) + pt(k+1,j,i) ) &
279	- ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
280	) / pt_surface
281	ENDDO
282	ENDDO
283	ENDDO
284
285	ELSE
286
287	WRITE( message_string, * ) 'no term for component "', &
288	wind_component,'"'
289	CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
290
291	ENDIF
292
293	ENDIF
294
295	END SUBROUTINE buoyancy_acc
296
297
298	!------------------------------------------------------------------------------!
299	! Call for grid point i,j
300	! ATTENTION: PGI-compiler creates SIGFPE if opt>1 is used! Therefore, opt=1 is
301	! forced by compiler-directive.
302	!------------------------------------------------------------------------------!
303	!pgi$r opt=1
304	SUBROUTINE buoyancy_ij( i, j, var, wind_component )
305
306	USE arrays_3d, &
307	ONLY: pt, pt_slope_ref, ref_state, tend
308
309	USE control_parameters, &
310	ONLY: atmos_ocean_sign, cos_alpha_surface, g, message_string, &
311	pt_surface, sin_alpha_surface, sloping_surface
312
313	USE indices, &
314	ONLY: nxlg, nxrg, nyng, nysg, nzb, nzb_s_inner, nzt
315
316	USE kinds
317
318	USE pegrid
319
320
321	IMPLICIT NONE
322
323	INTEGER(iwp) :: i !:
324	INTEGER(iwp) :: j !:
325	INTEGER(iwp) :: k !:
326	INTEGER(iwp) :: pr !:
327	INTEGER(iwp) :: wind_component !:
328
329	#if defined( __nopointer )
330	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !:
331	#else
332	REAL(wp), DIMENSION(:,:,:), POINTER :: var
333	#endif
334
335
336	IF ( .NOT. sloping_surface ) THEN
337	!
338	!-- Normal case: horizontal surface
339	DO k = nzb_s_inner(j,i)+1, nzt-1
340	tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5_wp * ( &
341	( var(k,j,i) - ref_state(k) ) / ref_state(k) + &
342	( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1) &
343	)
344	ENDDO
345
346	ELSE
347	!
348	!-- Buoyancy term for a surface with a slope in x-direction. The equations
349	!-- for both the u and w velocity-component contain proportionate terms.
350	!-- Temperature field at time t=0 serves as environmental temperature.
351	!-- Reference temperature (pt_surface) is the one at the lower left corner
352	!-- of the total domain.
353	IF ( wind_component == 1 ) THEN
354
355	DO k = nzb_s_inner(j,i)+1, nzt-1
356	tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface * &
357	0.5_wp * ( ( pt(k,j,i-1) + pt(k,j,i) ) &
358	- ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
359	) / pt_surface
360	ENDDO
361
362	ELSEIF ( wind_component == 3 ) THEN
363
364	DO k = nzb_s_inner(j,i)+1, nzt-1
365	tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface * &
366	0.5_wp * ( ( pt(k,j,i) + pt(k+1,j,i) ) &
367	- ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
368	) / pt_surface
369	ENDDO
370
371	ELSE
372
373	WRITE( message_string, * ) 'no term for component "', &
374	wind_component,'"'
375	CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
376
377	ENDIF
378
379	ENDIF
380
381	END SUBROUTINE buoyancy_ij
382
383	END MODULE buoyancy_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |