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

Last change on this file since 1171 was 1171, checked in by raasch, 8 years ago

New:
---

use_reference-case activated in accelerator version. (buoyancy, diffusion_e)
new option -e which defines the execution command to be used to run PALM,
compiler options for pgi/openacc added (palm_simple_run)
parameter sets for openACC benchmarks added (trunk/EXAMPLES/benchmark_acc)

Changed:


split of prognostic_equations deactivated (time_integration)

Errors:


bugfix: diss array is allocated with full size if accelerator boards are used (init_3d_model)

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