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

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

last commit documented

  • Property svn:keywords set to Id
File size: 15.7 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!
23!
24! Former revisions:
25! -----------------
26! $Id: buoyancy.f90 1154 2013-05-10 14:43:15Z 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             DO  i = i_left, i_right
228                DO  j = j_south, j_north
229                   DO  k = nzb_s_inner(j,i)+1, nzt-1
230                      tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * &
231                                                            (                  &
232                          ( var(k,j,i)   - hom(k,1,pr,0)   ) / var_reference + &
233                          ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / var_reference   &
234                                                            )
235                   ENDDO
236                ENDDO
237             ENDDO
238          ELSE
239             !$acc kernels present( nzb_s_inner, hom, tend, var )
240             !$acc loop
241             DO  i = i_left, i_right
242                DO  j = j_south, j_north
243                   !$acc loop independent vector(32)
244                   DO  k = nzb_s_inner(j,i)+1, nzt-1
245                      tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * &
246                                                            (                  &
247                          ( var(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
248                          ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
249                                                            )
250                   ENDDO
251                ENDDO
252             ENDDO
253             !$acc end kernels
254          ENDIF
255
256       ELSE
257!
258!--       Buoyancy term for a surface with a slope in x-direction. The equations
259!--       for both the u and w velocity-component contain proportionate terms.
260!--       Temperature field at time t=0 serves as environmental temperature.
261!--       Reference temperature (pt_surface) is the one at the lower left corner
262!--       of the total domain.
263          IF ( wind_component == 1 )  THEN
264
265             DO  i = nxlu, nxr
266                DO  j = nys, nyn
267                   DO  k = nzb_s_inner(j,i)+1, nzt-1
268                      tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *      &
269                           0.5 * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
270                                 - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
271                                 ) / pt_surface
272                   ENDDO
273                ENDDO
274             ENDDO
275
276          ELSEIF ( wind_component == 3 )  THEN
277
278             DO  i = nxl, nxr
279                DO  j = nys, nyn
280                   DO  k = nzb_s_inner(j,i)+1, nzt-1
281                      tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *      &
282                           0.5 * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
283                                 - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
284                                 ) / pt_surface
285                   ENDDO
286                ENDDO
287            ENDDO
288
289          ELSE
290
291             WRITE( message_string, * ) 'no term for component "',&
292                                       wind_component,'"'
293             CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
294
295          ENDIF
296
297       ENDIF
298
299    END SUBROUTINE buoyancy_acc
300
301
302!------------------------------------------------------------------------------!
303! Call for grid point i,j
304! ATTENTION: PGI-compiler creates SIGFPE if opt>1 is used! Therefore, opt=1 is
305!            forced by compiler-directive.
306!------------------------------------------------------------------------------!
307!pgi$r opt=1
308    SUBROUTINE buoyancy_ij( i, j, var, var_reference, wind_component, pr )
309
310       USE arrays_3d
311       USE control_parameters
312       USE indices
313       USE pegrid
314       USE statistics
315
316       IMPLICIT NONE
317
318       INTEGER ::  i, j, k, pr, wind_component
319       REAL    ::  var_reference
320#if defined( __nopointer )
321       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var
322#else
323       REAL, DIMENSION(:,:,:), POINTER ::  var
324#endif
325
326
327       IF ( .NOT. sloping_surface )  THEN
328!
329!--       Normal case: horizontal surface
330          IF ( use_reference )  THEN
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 * (   &
333                         ( var(k,j,i)   - hom(k,1,pr,0)   ) / var_reference + &
334                         ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / var_reference   &
335                                                                          )
336             ENDDO
337          ELSE
338             DO  k = nzb_s_inner(j,i)+1, nzt-1
339                 tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * (    &
340                          ( var(k,j,i)   - hom(k,1,pr,0)   ) / hom(k,1,pr,0) + &
341                          ( var(k+1,j,i) - hom(k+1,1,pr,0) ) / hom(k+1,1,pr,0) &
342                                                                          )
343             ENDDO
344          ENDIF
345
346       ELSE
347!
348!--       Buoyancy term for a surface with a slope in x-direction. The equations
349!--       for both the u and w velocity-component contain proportionate terms.
350!--       Temperature field at time t=0 serves as environmental temperature.
351!--       Reference temperature (pt_surface) is the one at the lower left corner
352!--       of the total domain.
353          IF ( wind_component == 1 )  THEN
354
355             DO  k = nzb_s_inner(j,i)+1, nzt-1
356                tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *            &
357                           0.5 * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
358                                 - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
359                                 ) / pt_surface
360             ENDDO
361
362          ELSEIF ( wind_component == 3 )  THEN
363
364             DO  k = nzb_s_inner(j,i)+1, nzt-1
365                tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *            &
366                           0.5 * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
367                                 - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
368                                 ) / pt_surface
369             ENDDO
370
371          ELSE
372
373             WRITE( message_string, * ) 'no term for component "',&
374                                       wind_component,'"'
375             CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
376
377          ENDIF
378
379       ENDIF
380
381    END SUBROUTINE buoyancy_ij
382
383
384    SUBROUTINE calc_mean_profile( var, pr )
385
386!------------------------------------------------------------------------------!
387! Description:
388! ------------
389! Calculate the horizontally averaged vertical temperature profile (pr=4 in case
390! of potential temperature and 44 in case of virtual potential temperature).
391!------------------------------------------------------------------------------!
392
393       USE control_parameters
394       USE indices
395       USE pegrid
396       USE statistics
397
398       IMPLICIT NONE
399
400       INTEGER ::  i, j, k, omp_get_thread_num, pr, tn
401#if defined( __nopointer )
402       REAL, DIMENSION(:,:,:) ::  var
403#else
404       REAL, DIMENSION(:,:,:), POINTER ::  var
405#endif
406
407!
408!--    Computation of the horizontally averaged profile of variable var, unless
409!--    already done by the relevant call from flow_statistics. The calculation
410!--    is done only for the first respective intermediate timestep in order to
411!--    spare communication time and to produce identical model results with jobs
412!--    which are calling flow_statistics at different time intervals.
413       IF ( .NOT. flow_statistics_called  .AND. &
414            intermediate_timestep_count == 1 )  THEN
415
416!
417!--       Horizontal average of variable var
418          tn           =   0  ! Default thread number in case of one thread
419          !$OMP PARALLEL PRIVATE( i, j, k, tn )
420!$        tn = omp_get_thread_num()
421          sums_l(:,pr,tn) = 0.0
422          !$OMP DO
423          DO  i = nxl, nxr
424             DO  j =  nys, nyn
425                DO  k = nzb_s_inner(j,i), nzt+1
426                   sums_l(k,pr,tn) = sums_l(k,pr,tn) + var(k,j,i)
427                ENDDO
428             ENDDO
429          ENDDO
430          !$OMP END PARALLEL
431
432          DO  i = 1, threads_per_task-1
433             sums_l(:,pr,0) = sums_l(:,pr,0) + sums_l(:,pr,i)
434          ENDDO
435
436#if defined( __parallel )
437
438          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
439          CALL MPI_ALLREDUCE( sums_l(nzb,pr,0), sums(nzb,pr), nzt+2-nzb, &
440                              MPI_REAL, MPI_SUM, comm2d, ierr )
441
442#else
443
444          sums(:,pr) = sums_l(:,pr,0)
445
446#endif
447
448          hom(:,1,pr,0) = sums(:,pr) / ngp_2dh_s_inner(:,0)
449
450       ENDIF
451
452    END SUBROUTINE calc_mean_profile
453
454 END MODULE buoyancy_mod
Note: See TracBrowser for help on using the repository browser.