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

Last change on this file since 1153 was 1153, checked in by raasch, 9 years ago

code adjustment of accelerator version for PGI 12.3 / CUDA 5.0

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