Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |