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

Last change on this file since 1374 was 1374, checked in by raasch, 10 years ago

bugfixes:
missing variables added to ONLY lists in USE statements (advec_s_bc, advec_s_pw, advec_s_up, advec_ws, buoyancy, diffusion_e, diffusion_s, fft_xy, flow_statistics, palm, prognostic_equations)
missing module kinds added (cuda_fft_interfaces)
dpk renamed dp (fft_xy)
missing dependency for check_open added (Makefile)
variables removed from acc-present-list (diffusion_e, diffusion_w, diffusivities, production_e, wall_fluxes)
syntax errors removed from openacc-branch (flow_statistics)
USE-statement for nopointer-case added (swap_timelevel)

  • Property svn:keywords set to Id
File size: 13.4 KB
RevLine 
[1]1 MODULE buoyancy_mod
2
[1036]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!
[1310]17! Copyright 1997-2014 Leibniz Universitaet Hannover
[1036]18!--------------------------------------------------------------------------------!
19!
[1328]20! Current revisions:
[1179]21! ------------------
[1374]22! missing variables added to ONLY list
[1354]23!
[1321]24! Former revisions:
25! -----------------
26! $Id: buoyancy.f90 1374 2014-04-25 12:55:07Z raasch $
27!
[1366]28! 1365 2014-04-22 15:03:56Z boeske
29! Calculation of reference state in subroutine calc_mean_profile moved to
30! subroutine time_integration,
31! subroutine calc_mean_profile moved to new file calc_mean_profile.f90
32!
[1354]33! 1353 2014-04-08 15:21:23Z heinze
34! REAL constants provided with KIND-attribute
35!
[1321]36! 1320 2014-03-20 08:40:49Z raasch
[1320]37! ONLY-attribute added to USE-statements,
38! kind-parameters added to all INTEGER and REAL declaration statements,
39! kinds are defined in new module kinds,
40! revision history before 2012 removed,
41! comment fields (!:) to be used for variable explanations added to
42! all variable declaration statements
[98]43!
[1258]44! 1257 2013-11-08 15:18:40Z raasch
45! vector length (32) removed from openacc clause
46!
[1242]47! 1241 2013-10-30 11:36:58Z heinze
48! Generalize calc_mean_profile for wider use: use additional steering
49! character loc
50!
[1182]51! 1179 2013-06-14 05:57:58Z raasch
52! steering of reference state revised (var_reference and pr removed from
53! parameter list), use_reference renamed use_single_reference_value
54!
[1172]55! 1171 2013-05-30 11:27:45Z raasch
56! openacc statements added to use_reference-case in accelerator version
57!
[1154]58! 1153 2013-05-10 14:33:08Z raasch
59! code adjustments of accelerator version required by PGI 12.3 / CUDA 5.0
60!
[1132]61! 1128 2013-04-12 06:19:32Z raasch
62! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
63! j_north
64!
[1037]65! 1036 2012-10-22 13:43:42Z raasch
66! code put under GPL (PALM 3.9)
67!
[1017]68! 1015 2012-09-27 09:23:24Z raasch
69! accelerator version (*_acc) added
70!
[1011]71! 1010 2012-09-20 07:59:54Z raasch
72! cpp switch __nopointer added for pointer free version
73!
[1]74! Revision 1.1  1997/08/29 08:56:48  raasch
75! Initial revision
76!
77!
78! Description:
79! ------------
80! Buoyancy term of the third component of the equation of motion.
81! WARNING: humidity is not regarded when using a sloping surface!
82!------------------------------------------------------------------------------!
83
84    PRIVATE
[1365]85    PUBLIC buoyancy, buoyancy_acc
[1]86
87    INTERFACE buoyancy
88       MODULE PROCEDURE buoyancy
89       MODULE PROCEDURE buoyancy_ij
90    END INTERFACE buoyancy
91
[1015]92    INTERFACE buoyancy_acc
93       MODULE PROCEDURE buoyancy_acc
94    END INTERFACE buoyancy_acc
95
[1]96 CONTAINS
97
98
99!------------------------------------------------------------------------------!
100! Call for all grid points
101!------------------------------------------------------------------------------!
[1179]102    SUBROUTINE buoyancy( var, wind_component )
[1]103
[1320]104       USE arrays_3d,                                                          &
105           ONLY:  pt, pt_slope_ref, ref_state, tend
106
107       USE control_parameters,                                                 &
108           ONLY:  atmos_ocean_sign, cos_alpha_surface, g, message_string,      &
109                  pt_surface, sin_alpha_surface, sloping_surface
110
111       USE indices,                                                            &
[1374]112           ONLY:  nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nzb,       &
113                  nzb_s_inner, nzt
[1320]114
115       USE kinds
116
[1]117       USE pegrid
118
[1320]119
[1]120       IMPLICIT NONE
121
[1320]122       INTEGER(iwp) ::  i              !:
123       INTEGER(iwp) ::  j              !:
124       INTEGER(iwp) ::  k              !:
125       INTEGER(iwp) ::  wind_component !:
126       
[1010]127#if defined( __nopointer )
[1320]128       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var !:
[1010]129#else
[1320]130       REAL(wp), DIMENSION(:,:,:), POINTER ::  var
[1010]131#endif
[1]132
133
134       IF ( .NOT. sloping_surface )  THEN
135!
136!--       Normal case: horizontal surface
[1179]137          DO  i = nxl, nxr
138             DO  j = nys, nyn
139                DO  k = nzb_s_inner(j,i)+1, nzt-1
[1353]140                   tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5_wp * &
141                          (                                                    &
142                             ( var(k,j,i)   - ref_state(k) )   / ref_state(k) + &
143                             ( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1) &
144                          )
[57]145                ENDDO
146             ENDDO
[1179]147          ENDDO
[1]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
[106]158             DO  i = nxlu, nxr
[1]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 *      &
[1353]162                           0.5_wp * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
163                                    - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
164                                    ) / pt_surface
[1]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 *      &
[1353]175                           0.5_wp * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
176                                    - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
177                                    ) / pt_surface
[1]178                   ENDDO
179                ENDDO
180            ENDDO
181
182          ELSE
[247]183             
[1320]184             WRITE( message_string, * ) 'no term for component "',             &
[1]185                                       wind_component,'"'
[247]186             CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
[1]187
188          ENDIF
189
190       ENDIF
191
192    END SUBROUTINE buoyancy
193
194
195!------------------------------------------------------------------------------!
[1015]196! Call for all grid points - accelerator version
197!------------------------------------------------------------------------------!
[1179]198    SUBROUTINE buoyancy_acc( var, wind_component )
[1015]199
[1320]200       USE arrays_3d,                                                          &
201           ONLY:  pt, pt_slope_ref, ref_state, tend
202
203       USE control_parameters,                                                 &
204           ONLY:  atmos_ocean_sign, cos_alpha_surface, g, message_string,      &
205                  pt_surface, sin_alpha_surface, sloping_surface
206
207       USE indices,                                                            &
[1374]208           ONLY:  i_left, i_right, j_north, j_south, nxl, nxlg, nxlu, nxr,     &
209                  nxrg, nyn, nyng, nys, nysg, nzb, nzb_s_inner, nzt
[1320]210
211       USE kinds
212
[1015]213       USE pegrid
214
[1320]215
[1015]216       IMPLICIT NONE
217
[1320]218       INTEGER(iwp) ::  i              !:
219       INTEGER(iwp) ::  j              !:
220       INTEGER(iwp) ::  k              !:
221       INTEGER(iwp) ::  wind_component !:
222       
[1015]223#if defined( __nopointer )
[1320]224       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var !:
[1015]225#else
[1320]226       REAL(wp), DIMENSION(:,:,:), POINTER ::  var
[1015]227#endif
228
229
230       IF ( .NOT. sloping_surface )  THEN
231!
232!--       Normal case: horizontal surface
[1179]233          !$acc kernels present( nzb_s_inner, ref_state, tend, var )
234          !$acc loop
235          DO  i = i_left, i_right
236             DO  j = j_south, j_north
[1257]237                !$acc loop independent vector
[1179]238                DO  k = nzb_s_inner(j,i)+1, nzt-1
[1353]239                   tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5_wp * &
240                          (                                                    &
241                             ( var(k,j,i)   - ref_state(k) )   / ref_state(k) + &
242                             ( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1) &
243                          )
[1015]244                ENDDO
245             ENDDO
[1179]246          ENDDO
247          !$acc end kernels
[1015]248
249       ELSE
250!
251!--       Buoyancy term for a surface with a slope in x-direction. The equations
252!--       for both the u and w velocity-component contain proportionate terms.
253!--       Temperature field at time t=0 serves as environmental temperature.
254!--       Reference temperature (pt_surface) is the one at the lower left corner
255!--       of the total domain.
256          IF ( wind_component == 1 )  THEN
257
258             DO  i = nxlu, nxr
259                DO  j = nys, nyn
260                   DO  k = nzb_s_inner(j,i)+1, nzt-1
261                      tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *      &
[1353]262                           0.5_wp * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
263                                    - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
264                                    ) / pt_surface
[1015]265                   ENDDO
266                ENDDO
267             ENDDO
268
269          ELSEIF ( wind_component == 3 )  THEN
270
271             DO  i = nxl, nxr
272                DO  j = nys, nyn
273                   DO  k = nzb_s_inner(j,i)+1, nzt-1
274                      tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *      &
[1353]275                           0.5_wp * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
276                                    - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
277                                    ) / pt_surface
[1015]278                   ENDDO
279                ENDDO
280            ENDDO
281
282          ELSE
283
[1320]284             WRITE( message_string, * ) 'no term for component "',             &
[1015]285                                       wind_component,'"'
286             CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
287
288          ENDIF
289
290       ENDIF
291
292    END SUBROUTINE buoyancy_acc
293
294
295!------------------------------------------------------------------------------!
[1]296! Call for grid point i,j
[515]297! ATTENTION: PGI-compiler creates SIGFPE if opt>1 is used! Therefore, opt=1 is
298!            forced by compiler-directive.
[1]299!------------------------------------------------------------------------------!
[515]300!pgi$r opt=1
[1179]301    SUBROUTINE buoyancy_ij( i, j, var, wind_component )
[1]302
[1320]303       USE arrays_3d,                                                          &
304           ONLY:  pt, pt_slope_ref, ref_state, tend
305
306       USE control_parameters,                                                 &
307           ONLY:  atmos_ocean_sign, cos_alpha_surface, g, message_string,      &
308                  pt_surface, sin_alpha_surface, sloping_surface
309
310       USE indices,                                                            &
[1374]311           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_s_inner, nzt
[1320]312
313       USE kinds
314
[1]315       USE pegrid
316
[1320]317
[1]318       IMPLICIT NONE
319
[1320]320       INTEGER(iwp) ::  i              !:
321       INTEGER(iwp) ::  j              !:
322       INTEGER(iwp) ::  k              !:
323       INTEGER(iwp) ::  pr             !:
324       INTEGER(iwp) ::  wind_component !:
325       
[1010]326#if defined( __nopointer )
[1320]327       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var !:
[1010]328#else
[1320]329       REAL(wp), DIMENSION(:,:,:), POINTER ::  var
[1010]330#endif
[1]331
332
333       IF ( .NOT. sloping_surface )  THEN
334!
335!--       Normal case: horizontal surface
[1179]336          DO  k = nzb_s_inner(j,i)+1, nzt-1
[1353]337              tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5_wp * (    &
[1320]338                        ( var(k,j,i)   - ref_state(k)   ) / ref_state(k)   +   &
339                        ( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1)     &
[1353]340                                                                          )
[1179]341          ENDDO
[1]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 *            &
[1353]354                           0.5_wp * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
355                                    - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
356                                    ) / pt_surface
[1]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 *            &
[1353]363                           0.5_wp * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
364                                    - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
365                                    ) / pt_surface
[1]366             ENDDO
367
368          ELSE
369
[1320]370             WRITE( message_string, * ) 'no term for component "',             &
[1]371                                       wind_component,'"'
[247]372             CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
[1]373
374          ENDIF
375
376       ENDIF
377
378    END SUBROUTINE buoyancy_ij
379
380 END MODULE buoyancy_mod
Note: See TracBrowser for help on using the repository browser.