Home

Context Navigation

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

Last change on this file since 1353 was 1353, checked in by heinze, 10 years ago
REAL constants provided with KIND-attribute
Property svn:keywords set to `Id`
File size: 16.6 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	! REAL constants provided with KIND-attribute
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: buoyancy.f90 1353 2014-04-08 15:21:23Z heinze $
27	!
28	! 1320 2014-03-20 08:40:49Z raasch
29	! ONLY-attribute added to USE-statements,
30	! kind-parameters added to all INTEGER and REAL declaration statements,
31	! kinds are defined in new module kinds,
32	! revision history before 2012 removed,
33	! comment fields (!:) to be used for variable explanations added to
34	! all variable declaration statements
35	!
36	! 1257 2013-11-08 15:18:40Z raasch
37	! vector length (32) removed from openacc clause
38	!
39	! 1241 2013-10-30 11:36:58Z heinze
40	! Generalize calc_mean_profile for wider use: use additional steering
41	! character loc
42	!
43	! 1179 2013-06-14 05:57:58Z raasch
44	! steering of reference state revised (var_reference and pr removed from
45	! parameter list), use_reference renamed use_single_reference_value
46	!
47	! 1171 2013-05-30 11:27:45Z raasch
48	! openacc statements added to use_reference-case in accelerator version
49	!
50	! 1153 2013-05-10 14:33:08Z raasch
51	! code adjustments of accelerator version required by PGI 12.3 / CUDA 5.0
52	!
53	! 1128 2013-04-12 06:19:32Z raasch
54	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
55	! j_north
56	!
57	! 1036 2012-10-22 13:43:42Z raasch
58	! code put under GPL (PALM 3.9)
59	!
60	! 1015 2012-09-27 09:23:24Z raasch
61	! accelerator version (*_acc) added
62	!
63	! 1010 2012-09-20 07:59:54Z raasch
64	! cpp switch __nopointer added for pointer free version
65	!
66	! Revision 1.1 1997/08/29 08:56:48 raasch
67	! Initial revision
68	!
69	!
70	! Description:
71	! ------------
72	! Buoyancy term of the third component of the equation of motion.
73	! WARNING: humidity is not regarded when using a sloping surface!
74	!------------------------------------------------------------------------------!
75
76	PRIVATE
77	PUBLIC buoyancy, buoyancy_acc, calc_mean_profile
78
79	INTERFACE buoyancy
80	MODULE PROCEDURE buoyancy
81	MODULE PROCEDURE buoyancy_ij
82	END INTERFACE buoyancy
83
84	INTERFACE buoyancy_acc
85	MODULE PROCEDURE buoyancy_acc
86	END INTERFACE buoyancy_acc
87
88	INTERFACE calc_mean_profile
89	MODULE PROCEDURE calc_mean_profile
90	END INTERFACE calc_mean_profile
91
92	CONTAINS
93
94
95	!------------------------------------------------------------------------------!
96	! Call for all grid points
97	!------------------------------------------------------------------------------!
98	SUBROUTINE buoyancy( var, wind_component )
99
100	USE arrays_3d, &
101	ONLY: pt, pt_slope_ref, ref_state, tend
102
103	USE control_parameters, &
104	ONLY: atmos_ocean_sign, cos_alpha_surface, g, message_string, &
105	pt_surface, sin_alpha_surface, sloping_surface
106
107	USE indices, &
108	ONLY: nxl, nxlu, nxr, nyn, nys, nzb_s_inner, nzt
109
110	USE kinds
111
112	USE pegrid
113
114
115	IMPLICIT NONE
116
117	INTEGER(iwp) :: i !:
118	INTEGER(iwp) :: j !:
119	INTEGER(iwp) :: k !:
120	INTEGER(iwp) :: wind_component !:
121
122	#if defined( __nopointer )
123	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !:
124	#else
125	REAL(wp), DIMENSION(:,:,:), POINTER :: var
126	#endif
127
128
129	IF ( .NOT. sloping_surface ) THEN
130	!
131	!-- Normal case: horizontal surface
132	DO i = nxl, nxr
133	DO j = nys, nyn
134	DO k = nzb_s_inner(j,i)+1, nzt-1
135	tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5_wp * &
136	( &
137	( var(k,j,i) - ref_state(k) ) / ref_state(k) + &
138	( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1) &
139	)
140	ENDDO
141	ENDDO
142	ENDDO
143
144	ELSE
145	!
146	!-- Buoyancy term for a surface with a slope in x-direction. The equations
147	!-- for both the u and w velocity-component contain proportionate terms.
148	!-- Temperature field at time t=0 serves as environmental temperature.
149	!-- Reference temperature (pt_surface) is the one at the lower left corner
150	!-- of the total domain.
151	IF ( wind_component == 1 ) THEN
152
153	DO i = nxlu, nxr
154	DO j = nys, nyn
155	DO k = nzb_s_inner(j,i)+1, nzt-1
156	tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface * &
157	0.5_wp * ( ( pt(k,j,i-1) + pt(k,j,i) ) &
158	- ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
159	) / pt_surface
160	ENDDO
161	ENDDO
162	ENDDO
163
164	ELSEIF ( wind_component == 3 ) THEN
165
166	DO i = nxl, nxr
167	DO j = nys, nyn
168	DO k = nzb_s_inner(j,i)+1, nzt-1
169	tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface * &
170	0.5_wp * ( ( pt(k,j,i) + pt(k+1,j,i) ) &
171	- ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
172	) / pt_surface
173	ENDDO
174	ENDDO
175	ENDDO
176
177	ELSE
178
179	WRITE( message_string, * ) 'no term for component "', &
180	wind_component,'"'
181	CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
182
183	ENDIF
184
185	ENDIF
186
187	END SUBROUTINE buoyancy
188
189
190	!------------------------------------------------------------------------------!
191	! Call for all grid points - accelerator version
192	!------------------------------------------------------------------------------!
193	SUBROUTINE buoyancy_acc( var, wind_component )
194
195	USE arrays_3d, &
196	ONLY: pt, pt_slope_ref, ref_state, tend
197
198	USE control_parameters, &
199	ONLY: atmos_ocean_sign, cos_alpha_surface, g, message_string, &
200	pt_surface, sin_alpha_surface, sloping_surface
201
202	USE indices, &
203	ONLY: i_left, i_right, j_north, j_south, nxl, nxlu, nxr, nyn, nys, &
204	nzb_s_inner, nzt
205
206	USE kinds
207
208	USE pegrid
209
210
211	IMPLICIT NONE
212
213	INTEGER(iwp) :: i !:
214	INTEGER(iwp) :: j !:
215	INTEGER(iwp) :: k !:
216	INTEGER(iwp) :: wind_component !:
217
218	#if defined( __nopointer )
219	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !:
220	#else
221	REAL(wp), DIMENSION(:,:,:), POINTER :: var
222	#endif
223
224
225	IF ( .NOT. sloping_surface ) THEN
226	!
227	!-- Normal case: horizontal surface
228	!$acc kernels present( nzb_s_inner, ref_state, tend, var )
229	!$acc loop
230	DO i = i_left, i_right
231	DO j = j_south, j_north
232	!$acc loop independent vector
233	DO k = nzb_s_inner(j,i)+1, nzt-1
234	tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5_wp * &
235	( &
236	( var(k,j,i) - ref_state(k) ) / ref_state(k) + &
237	( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1) &
238	)
239	ENDDO
240	ENDDO
241	ENDDO
242	!$acc end kernels
243
244	ELSE
245	!
246	!-- Buoyancy term for a surface with a slope in x-direction. The equations
247	!-- for both the u and w velocity-component contain proportionate terms.
248	!-- Temperature field at time t=0 serves as environmental temperature.
249	!-- Reference temperature (pt_surface) is the one at the lower left corner
250	!-- of the total domain.
251	IF ( wind_component == 1 ) THEN
252
253	DO i = nxlu, nxr
254	DO j = nys, nyn
255	DO k = nzb_s_inner(j,i)+1, nzt-1
256	tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface * &
257	0.5_wp * ( ( pt(k,j,i-1) + pt(k,j,i) ) &
258	- ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
259	) / pt_surface
260	ENDDO
261	ENDDO
262	ENDDO
263
264	ELSEIF ( wind_component == 3 ) THEN
265
266	DO i = nxl, nxr
267	DO j = nys, nyn
268	DO k = nzb_s_inner(j,i)+1, nzt-1
269	tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface * &
270	0.5_wp * ( ( pt(k,j,i) + pt(k+1,j,i) ) &
271	- ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
272	) / pt_surface
273	ENDDO
274	ENDDO
275	ENDDO
276
277	ELSE
278
279	WRITE( message_string, * ) 'no term for component "', &
280	wind_component,'"'
281	CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
282
283	ENDIF
284
285	ENDIF
286
287	END SUBROUTINE buoyancy_acc
288
289
290	!------------------------------------------------------------------------------!
291	! Call for grid point i,j
292	! ATTENTION: PGI-compiler creates SIGFPE if opt>1 is used! Therefore, opt=1 is
293	! forced by compiler-directive.
294	!------------------------------------------------------------------------------!
295	!pgi$r opt=1
296	SUBROUTINE buoyancy_ij( i, j, var, wind_component )
297
298	USE arrays_3d, &
299	ONLY: pt, pt_slope_ref, ref_state, tend
300
301	USE control_parameters, &
302	ONLY: atmos_ocean_sign, cos_alpha_surface, g, message_string, &
303	pt_surface, sin_alpha_surface, sloping_surface
304
305	USE indices, &
306	ONLY: nzb_s_inner, nzt
307
308	USE kinds
309
310	USE pegrid
311
312
313	IMPLICIT NONE
314
315	INTEGER(iwp) :: i !:
316	INTEGER(iwp) :: j !:
317	INTEGER(iwp) :: k !:
318	INTEGER(iwp) :: pr !:
319	INTEGER(iwp) :: wind_component !:
320
321	#if defined( __nopointer )
322	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !:
323	#else
324	REAL(wp), DIMENSION(:,:,:), POINTER :: var
325	#endif
326
327
328	IF ( .NOT. sloping_surface ) THEN
329	!
330	!-- Normal case: horizontal surface
331	DO k = nzb_s_inner(j,i)+1, nzt-1
332	tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5_wp * ( &
333	( var(k,j,i) - ref_state(k) ) / ref_state(k) + &
334	( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1) &
335	)
336	ENDDO
337
338	ELSE
339	!
340	!-- Buoyancy term for a surface with a slope in x-direction. The equations
341	!-- for both the u and w velocity-component contain proportionate terms.
342	!-- Temperature field at time t=0 serves as environmental temperature.
343	!-- Reference temperature (pt_surface) is the one at the lower left corner
344	!-- of the total domain.
345	IF ( wind_component == 1 ) THEN
346
347	DO k = nzb_s_inner(j,i)+1, nzt-1
348	tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface * &
349	0.5_wp * ( ( pt(k,j,i-1) + pt(k,j,i) ) &
350	- ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
351	) / pt_surface
352	ENDDO
353
354	ELSEIF ( wind_component == 3 ) THEN
355
356	DO k = nzb_s_inner(j,i)+1, nzt-1
357	tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface * &
358	0.5_wp * ( ( pt(k,j,i) + pt(k+1,j,i) ) &
359	- ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
360	) / pt_surface
361	ENDDO
362
363	ELSE
364
365	WRITE( message_string, * ) 'no term for component "', &
366	wind_component,'"'
367	CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
368
369	ENDIF
370
371	ENDIF
372
373	END SUBROUTINE buoyancy_ij
374
375
376	SUBROUTINE calc_mean_profile( var, pr, loc )
377
378	!------------------------------------------------------------------------------!
379	! Description:
380	! ------------
381	! Calculate the horizontally averaged vertical temperature profile (pr=4 in case
382	! of potential temperature, 44 in case of virtual potential temperature, and 64
383	! in case of density (ocean runs)).
384	!------------------------------------------------------------------------------!
385
386	USE arrays_3d, &
387	ONLY: ref_state
388
389	USE control_parameters, &
390	ONLY: intermediate_timestep_count, message_string
391
392	USE indices, &
393	ONLY: ngp_2dh_s_inner, nxl, nxr, nyn, nys, nzb, nzb_s_inner, nzt
394
395	USE kinds
396
397	USE pegrid
398
399	USE statistics, &
400	ONLY: flow_statistics_called, hom, sums, sums_l
401
402
403	IMPLICIT NONE
404
405	CHARACTER (LEN=*) :: loc !:
406
407	INTEGER(iwp) :: i !:
408	INTEGER(iwp) :: j !:
409	INTEGER(iwp) :: k !:
410	INTEGER(iwp) :: pr !:
411	INTEGER(iwp) :: omp_get_thread_num !:
412	INTEGER(iwp) :: tn !:
413
414	#if defined( __nopointer )
415	REAL(wp), DIMENSION(:,:,:) :: var !:
416	#else
417	REAL(wp), DIMENSION(:,:,:), POINTER :: var
418	#endif
419
420	!
421	!-- Computation of the horizontally averaged profile of variable var, unless
422	!-- already done by the relevant call from flow_statistics. The calculation
423	!-- is done only for the first respective intermediate timestep in order to
424	!-- spare communication time and to produce identical model results with jobs
425	!-- which are calling flow_statistics at different time intervals.
426	IF ( .NOT. flow_statistics_called .AND. &
427	intermediate_timestep_count == 1 ) THEN
428
429	!
430	!-- Horizontal average of variable var
431	tn = 0 ! Default thread number in case of one thread
432	!$OMP PARALLEL PRIVATE( i, j, k, tn )
433	!$ tn = omp_get_thread_num()
434	sums_l(:,pr,tn) = 0.0_wp
435	!$OMP DO
436	DO i = nxl, nxr
437	DO j = nys, nyn
438	DO k = nzb_s_inner(j,i), nzt+1
439	sums_l(k,pr,tn) = sums_l(k,pr,tn) + var(k,j,i)
440	ENDDO
441	ENDDO
442	ENDDO
443	!$OMP END PARALLEL
444
445	DO i = 1, threads_per_task-1
446	sums_l(:,pr,0) = sums_l(:,pr,0) + sums_l(:,pr,i)
447	ENDDO
448
449	#if defined( __parallel )
450
451	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
452	CALL MPI_ALLREDUCE( sums_l(nzb,pr,0), sums(nzb,pr), nzt+2-nzb, &
453	MPI_REAL, MPI_SUM, comm2d, ierr )
454
455	#else
456
457	sums(:,pr) = sums_l(:,pr,0)
458
459	#endif
460
461	hom(:,1,pr,0) = sums(:,pr) / ngp_2dh_s_inner(:,0)
462
463	ENDIF
464
465	SELECT CASE ( loc )
466
467	CASE ( 'time_int' )
468
469	ref_state(:) = hom(:,1,pr,0) ! this is used in the buoyancy term
470
471
472	CASE ( 'nudging' )
473	!nothing to be done
474
475
476	CASE DEFAULT
477	message_string = 'unknown location "' // loc // '"'
478	CALL message( 'calc_mean_profile', 'PA0379', 1, 2, 0, 6, 0 )
479
480	END SELECT
481
482
483
484	END SUBROUTINE calc_mean_profile
485
486	END MODULE buoyancy_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |