Home

Context Navigation

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

Last change on this file since 1329 was 1328, checked in by maronga, 10 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 16.4 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 1328 2014-03-21 11:00:33Z raasch $
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 * ( &
136	( var(k,j,i) - ref_state(k) ) / ref_state(k) + &
137	( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1) &
138	)
139	ENDDO
140	ENDDO
141	ENDDO
142
143	ELSE
144	!
145	!-- Buoyancy term for a surface with a slope in x-direction. The equations
146	!-- for both the u and w velocity-component contain proportionate terms.
147	!-- Temperature field at time t=0 serves as environmental temperature.
148	!-- Reference temperature (pt_surface) is the one at the lower left corner
149	!-- of the total domain.
150	IF ( wind_component == 1 ) THEN
151
152	DO i = nxlu, nxr
153	DO j = nys, nyn
154	DO k = nzb_s_inner(j,i)+1, nzt-1
155	tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface * &
156	0.5 * ( ( pt(k,j,i-1) + pt(k,j,i) ) &
157	- ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
158	) / pt_surface
159	ENDDO
160	ENDDO
161	ENDDO
162
163	ELSEIF ( wind_component == 3 ) THEN
164
165	DO i = nxl, nxr
166	DO j = nys, nyn
167	DO k = nzb_s_inner(j,i)+1, nzt-1
168	tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface * &
169	0.5 * ( ( pt(k,j,i) + pt(k+1,j,i) ) &
170	- ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
171	) / pt_surface
172	ENDDO
173	ENDDO
174	ENDDO
175
176	ELSE
177
178	WRITE( message_string, * ) 'no term for component "', &
179	wind_component,'"'
180	CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
181
182	ENDIF
183
184	ENDIF
185
186	END SUBROUTINE buoyancy
187
188
189	!------------------------------------------------------------------------------!
190	! Call for all grid points - accelerator version
191	!------------------------------------------------------------------------------!
192	SUBROUTINE buoyancy_acc( var, wind_component )
193
194	USE arrays_3d, &
195	ONLY: pt, pt_slope_ref, ref_state, tend
196
197	USE control_parameters, &
198	ONLY: atmos_ocean_sign, cos_alpha_surface, g, message_string, &
199	pt_surface, sin_alpha_surface, sloping_surface
200
201	USE indices, &
202	ONLY: i_left, i_right, j_north, j_south, nxl, nxlu, nxr, nyn, nys, &
203	nzb_s_inner, nzt
204
205	USE kinds
206
207	USE pegrid
208
209
210	IMPLICIT NONE
211
212	INTEGER(iwp) :: i !:
213	INTEGER(iwp) :: j !:
214	INTEGER(iwp) :: k !:
215	INTEGER(iwp) :: wind_component !:
216
217	#if defined( __nopointer )
218	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !:
219	#else
220	REAL(wp), DIMENSION(:,:,:), POINTER :: var
221	#endif
222
223
224	IF ( .NOT. sloping_surface ) THEN
225	!
226	!-- Normal case: horizontal surface
227	!$acc kernels present( nzb_s_inner, ref_state, tend, var )
228	!$acc loop
229	DO i = i_left, i_right
230	DO j = j_south, j_north
231	!$acc loop independent vector
232	DO k = nzb_s_inner(j,i)+1, nzt-1
233	tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * ( &
234	( var(k,j,i) - ref_state(k) ) / ref_state(k) + &
235	( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1) &
236	)
237	ENDDO
238	ENDDO
239	ENDDO
240	!$acc end kernels
241
242	ELSE
243	!
244	!-- Buoyancy term for a surface with a slope in x-direction. The equations
245	!-- for both the u and w velocity-component contain proportionate terms.
246	!-- Temperature field at time t=0 serves as environmental temperature.
247	!-- Reference temperature (pt_surface) is the one at the lower left corner
248	!-- of the total domain.
249	IF ( wind_component == 1 ) THEN
250
251	DO i = nxlu, nxr
252	DO j = nys, nyn
253	DO k = nzb_s_inner(j,i)+1, nzt-1
254	tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface * &
255	0.5 * ( ( pt(k,j,i-1) + pt(k,j,i) ) &
256	- ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
257	) / pt_surface
258	ENDDO
259	ENDDO
260	ENDDO
261
262	ELSEIF ( wind_component == 3 ) THEN
263
264	DO i = nxl, nxr
265	DO j = nys, nyn
266	DO k = nzb_s_inner(j,i)+1, nzt-1
267	tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface * &
268	0.5 * ( ( pt(k,j,i) + pt(k+1,j,i) ) &
269	- ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
270	) / pt_surface
271	ENDDO
272	ENDDO
273	ENDDO
274
275	ELSE
276
277	WRITE( message_string, * ) 'no term for component "', &
278	wind_component,'"'
279	CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
280
281	ENDIF
282
283	ENDIF
284
285	END SUBROUTINE buoyancy_acc
286
287
288	!------------------------------------------------------------------------------!
289	! Call for grid point i,j
290	! ATTENTION: PGI-compiler creates SIGFPE if opt>1 is used! Therefore, opt=1 is
291	! forced by compiler-directive.
292	!------------------------------------------------------------------------------!
293	!pgi$r opt=1
294	SUBROUTINE buoyancy_ij( i, j, var, wind_component )
295
296	USE arrays_3d, &
297	ONLY: pt, pt_slope_ref, ref_state, tend
298
299	USE control_parameters, &
300	ONLY: atmos_ocean_sign, cos_alpha_surface, g, message_string, &
301	pt_surface, sin_alpha_surface, sloping_surface
302
303	USE indices, &
304	ONLY: nzb_s_inner, nzt
305
306	USE kinds
307
308	USE pegrid
309
310
311	IMPLICIT NONE
312
313	INTEGER(iwp) :: i !:
314	INTEGER(iwp) :: j !:
315	INTEGER(iwp) :: k !:
316	INTEGER(iwp) :: pr !:
317	INTEGER(iwp) :: wind_component !:
318
319	#if defined( __nopointer )
320	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !:
321	#else
322	REAL(wp), DIMENSION(:,:,:), POINTER :: var
323	#endif
324
325
326	IF ( .NOT. sloping_surface ) THEN
327	!
328	!-- Normal case: horizontal surface
329	DO k = nzb_s_inner(j,i)+1, nzt-1
330	tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * ( &
331	( var(k,j,i) - ref_state(k) ) / ref_state(k) + &
332	( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1) &
333	)
334	ENDDO
335
336	ELSE
337	!
338	!-- Buoyancy term for a surface with a slope in x-direction. The equations
339	!-- for both the u and w velocity-component contain proportionate terms.
340	!-- Temperature field at time t=0 serves as environmental temperature.
341	!-- Reference temperature (pt_surface) is the one at the lower left corner
342	!-- of the total domain.
343	IF ( wind_component == 1 ) THEN
344
345	DO k = nzb_s_inner(j,i)+1, nzt-1
346	tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface * &
347	0.5 * ( ( pt(k,j,i-1) + pt(k,j,i) ) &
348	- ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
349	) / pt_surface
350	ENDDO
351
352	ELSEIF ( wind_component == 3 ) THEN
353
354	DO k = nzb_s_inner(j,i)+1, nzt-1
355	tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface * &
356	0.5 * ( ( pt(k,j,i) + pt(k+1,j,i) ) &
357	- ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
358	) / pt_surface
359	ENDDO
360
361	ELSE
362
363	WRITE( message_string, * ) 'no term for component "', &
364	wind_component,'"'
365	CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
366
367	ENDIF
368
369	ENDIF
370
371	END SUBROUTINE buoyancy_ij
372
373
374	SUBROUTINE calc_mean_profile( var, pr, loc )
375
376	!------------------------------------------------------------------------------!
377	! Description:
378	! ------------
379	! Calculate the horizontally averaged vertical temperature profile (pr=4 in case
380	! of potential temperature, 44 in case of virtual potential temperature, and 64
381	! in case of density (ocean runs)).
382	!------------------------------------------------------------------------------!
383
384	USE arrays_3d, &
385	ONLY: ref_state
386
387	USE control_parameters, &
388	ONLY: intermediate_timestep_count, message_string
389
390	USE indices, &
391	ONLY: ngp_2dh_s_inner, nxl, nxr, nyn, nys, nzb, nzb_s_inner, nzt
392
393	USE kinds
394
395	USE pegrid
396
397	USE statistics, &
398	ONLY: flow_statistics_called, hom, sums, sums_l
399
400
401	IMPLICIT NONE
402
403	CHARACTER (LEN=*) :: loc !:
404
405	INTEGER(iwp) :: i !:
406	INTEGER(iwp) :: j !:
407	INTEGER(iwp) :: k !:
408	INTEGER(iwp) :: pr !:
409	INTEGER(iwp) :: omp_get_thread_num !:
410	INTEGER(iwp) :: tn !:
411
412	#if defined( __nopointer )
413	REAL(wp), DIMENSION(:,:,:) :: var !:
414	#else
415	REAL(wp), DIMENSION(:,:,:), POINTER :: var
416	#endif
417
418	!
419	!-- Computation of the horizontally averaged profile of variable var, unless
420	!-- already done by the relevant call from flow_statistics. The calculation
421	!-- is done only for the first respective intermediate timestep in order to
422	!-- spare communication time and to produce identical model results with jobs
423	!-- which are calling flow_statistics at different time intervals.
424	IF ( .NOT. flow_statistics_called .AND. &
425	intermediate_timestep_count == 1 ) THEN
426
427	!
428	!-- Horizontal average of variable var
429	tn = 0 ! Default thread number in case of one thread
430	!$OMP PARALLEL PRIVATE( i, j, k, tn )
431	!$ tn = omp_get_thread_num()
432	sums_l(:,pr,tn) = 0.0
433	!$OMP DO
434	DO i = nxl, nxr
435	DO j = nys, nyn
436	DO k = nzb_s_inner(j,i), nzt+1
437	sums_l(k,pr,tn) = sums_l(k,pr,tn) + var(k,j,i)
438	ENDDO
439	ENDDO
440	ENDDO
441	!$OMP END PARALLEL
442
443	DO i = 1, threads_per_task-1
444	sums_l(:,pr,0) = sums_l(:,pr,0) + sums_l(:,pr,i)
445	ENDDO
446
447	#if defined( __parallel )
448
449	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
450	CALL MPI_ALLREDUCE( sums_l(nzb,pr,0), sums(nzb,pr), nzt+2-nzb, &
451	MPI_REAL, MPI_SUM, comm2d, ierr )
452
453	#else
454
455	sums(:,pr) = sums_l(:,pr,0)
456
457	#endif
458
459	hom(:,1,pr,0) = sums(:,pr) / ngp_2dh_s_inner(:,0)
460
461	ENDIF
462
463	SELECT CASE ( loc )
464
465	CASE ( 'time_int' )
466
467	ref_state(:) = hom(:,1,pr,0) ! this is used in the buoyancy term
468
469
470	CASE ( 'nudging' )
471	!nothing to be done
472
473
474	CASE DEFAULT
475	message_string = 'unknown location "' // loc // '"'
476	CALL message( 'calc_mean_profile', 'PA0379', 1, 2, 0, 6, 0 )
477
478	END SELECT
479
480
481
482	END SUBROUTINE calc_mean_profile
483
484	END MODULE buoyancy_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |