Home

Context Navigation

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

Last change on this file since 1927 was 1874, checked in by maronga, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 13.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |