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

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

asynchronous transfer of ghost point data for acc-optimized version

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