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

Last change on this file since 1366 was 1366, checked in by boeske, 7 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 13.3 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-2014 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: buoyancy.f90 1366 2014-04-22 15:06:33Z boeske $
27!
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!
33! 1353 2014-04-08 15:21:23Z heinze
34! REAL constants provided with KIND-attribute
35!
36! 1320 2014-03-20 08:40:49Z raasch
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
43!
44! 1257 2013-11-08 15:18:40Z raasch
45! vector length (32) removed from openacc clause
46!
47! 1241 2013-10-30 11:36:58Z heinze
48! Generalize calc_mean_profile for wider use: use additional steering
49! character loc
50!
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!
55! 1171 2013-05-30 11:27:45Z raasch
56! openacc statements added to use_reference-case in accelerator version
57!
58! 1153 2013-05-10 14:33:08Z raasch
59! code adjustments of accelerator version required by PGI 12.3 / CUDA 5.0
60!
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!
65! 1036 2012-10-22 13:43:42Z raasch
66! code put under GPL (PALM 3.9)
67!
68! 1015 2012-09-27 09:23:24Z raasch
69! accelerator version (*_acc) added
70!
71! 1010 2012-09-20 07:59:54Z raasch
72! cpp switch __nopointer added for pointer free version
73!
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
85    PUBLIC buoyancy, buoyancy_acc
86
87    INTERFACE buoyancy
88       MODULE PROCEDURE buoyancy
89       MODULE PROCEDURE buoyancy_ij
90    END INTERFACE buoyancy
91
92    INTERFACE buoyancy_acc
93       MODULE PROCEDURE buoyancy_acc
94    END INTERFACE buoyancy_acc
95
96 CONTAINS
97
98
99!------------------------------------------------------------------------------!
100! Call for all grid points
101!------------------------------------------------------------------------------!
102    SUBROUTINE buoyancy( var, wind_component )
103
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,                                                            &
112           ONLY:  nxl, nxlu, nxr, nyn, nys, nzb_s_inner, nzt
113
114       USE kinds
115
116       USE pegrid
117
118
119       IMPLICIT NONE
120
121       INTEGER(iwp) ::  i              !:
122       INTEGER(iwp) ::  j              !:
123       INTEGER(iwp) ::  k              !:
124       INTEGER(iwp) ::  wind_component !:
125       
126#if defined( __nopointer )
127       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var !:
128#else
129       REAL(wp), DIMENSION(:,:,:), POINTER ::  var
130#endif
131
132
133       IF ( .NOT. sloping_surface )  THEN
134!
135!--       Normal case: horizontal surface
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_wp * &
140                          (                                                    &
141                             ( var(k,j,i)   - ref_state(k) )   / ref_state(k) + &
142                             ( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1) &
143                          )
144                ENDDO
145             ENDDO
146          ENDDO
147
148       ELSE
149!
150!--       Buoyancy term for a surface with a slope in x-direction. The equations
151!--       for both the u and w velocity-component contain proportionate terms.
152!--       Temperature field at time t=0 serves as environmental temperature.
153!--       Reference temperature (pt_surface) is the one at the lower left corner
154!--       of the total domain.
155          IF ( wind_component == 1 )  THEN
156
157             DO  i = nxlu, nxr
158                DO  j = nys, nyn
159                   DO  k = nzb_s_inner(j,i)+1, nzt-1
160                      tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *      &
161                           0.5_wp * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
162                                    - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
163                                    ) / pt_surface
164                   ENDDO
165                ENDDO
166             ENDDO
167
168          ELSEIF ( wind_component == 3 )  THEN
169
170             DO  i = nxl, nxr
171                DO  j = nys, nyn
172                   DO  k = nzb_s_inner(j,i)+1, nzt-1
173                      tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *      &
174                           0.5_wp * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
175                                    - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
176                                    ) / pt_surface
177                   ENDDO
178                ENDDO
179            ENDDO
180
181          ELSE
182             
183             WRITE( message_string, * ) 'no term for component "',             &
184                                       wind_component,'"'
185             CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
186
187          ENDIF
188
189       ENDIF
190
191    END SUBROUTINE buoyancy
192
193
194!------------------------------------------------------------------------------!
195! Call for all grid points - accelerator version
196!------------------------------------------------------------------------------!
197    SUBROUTINE buoyancy_acc( var, wind_component )
198
199       USE arrays_3d,                                                          &
200           ONLY:  pt, pt_slope_ref, ref_state, tend
201
202       USE control_parameters,                                                 &
203           ONLY:  atmos_ocean_sign, cos_alpha_surface, g, message_string,      &
204                  pt_surface, sin_alpha_surface, sloping_surface
205
206       USE indices,                                                            &
207           ONLY:  i_left, i_right, j_north, j_south, nxl, nxlu, nxr, nyn, nys, &
208                  nzb_s_inner, nzt
209
210       USE kinds
211
212       USE pegrid
213
214
215       IMPLICIT NONE
216
217       INTEGER(iwp) ::  i              !:
218       INTEGER(iwp) ::  j              !:
219       INTEGER(iwp) ::  k              !:
220       INTEGER(iwp) ::  wind_component !:
221       
222#if defined( __nopointer )
223       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var !:
224#else
225       REAL(wp), DIMENSION(:,:,:), POINTER ::  var
226#endif
227
228
229       IF ( .NOT. sloping_surface )  THEN
230!
231!--       Normal case: horizontal surface
232          !$acc kernels present( nzb_s_inner, ref_state, tend, var )
233          !$acc loop
234          DO  i = i_left, i_right
235             DO  j = j_south, j_north
236                !$acc loop independent vector
237                DO  k = nzb_s_inner(j,i)+1, nzt-1
238                   tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5_wp * &
239                          (                                                    &
240                             ( var(k,j,i)   - ref_state(k) )   / ref_state(k) + &
241                             ( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1) &
242                          )
243                ENDDO
244             ENDDO
245          ENDDO
246          !$acc end kernels
247
248       ELSE
249!
250!--       Buoyancy term for a surface with a slope in x-direction. The equations
251!--       for both the u and w velocity-component contain proportionate terms.
252!--       Temperature field at time t=0 serves as environmental temperature.
253!--       Reference temperature (pt_surface) is the one at the lower left corner
254!--       of the total domain.
255          IF ( wind_component == 1 )  THEN
256
257             DO  i = nxlu, nxr
258                DO  j = nys, nyn
259                   DO  k = nzb_s_inner(j,i)+1, nzt-1
260                      tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *      &
261                           0.5_wp * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
262                                    - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
263                                    ) / pt_surface
264                   ENDDO
265                ENDDO
266             ENDDO
267
268          ELSEIF ( wind_component == 3 )  THEN
269
270             DO  i = nxl, nxr
271                DO  j = nys, nyn
272                   DO  k = nzb_s_inner(j,i)+1, nzt-1
273                      tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *      &
274                           0.5_wp * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
275                                    - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
276                                    ) / pt_surface
277                   ENDDO
278                ENDDO
279            ENDDO
280
281          ELSE
282
283             WRITE( message_string, * ) 'no term for component "',             &
284                                       wind_component,'"'
285             CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
286
287          ENDIF
288
289       ENDIF
290
291    END SUBROUTINE buoyancy_acc
292
293
294!------------------------------------------------------------------------------!
295! Call for grid point i,j
296! ATTENTION: PGI-compiler creates SIGFPE if opt>1 is used! Therefore, opt=1 is
297!            forced by compiler-directive.
298!------------------------------------------------------------------------------!
299!pgi$r opt=1
300    SUBROUTINE buoyancy_ij( i, j, var, wind_component )
301
302       USE arrays_3d,                                                          &
303           ONLY:  pt, pt_slope_ref, ref_state, tend
304
305       USE control_parameters,                                                 &
306           ONLY:  atmos_ocean_sign, cos_alpha_surface, g, message_string,      &
307                  pt_surface, sin_alpha_surface, sloping_surface
308
309       USE indices,                                                            &
310           ONLY:  nzb_s_inner, nzt
311
312       USE kinds
313
314       USE pegrid
315
316
317       IMPLICIT NONE
318
319       INTEGER(iwp) ::  i              !:
320       INTEGER(iwp) ::  j              !:
321       INTEGER(iwp) ::  k              !:
322       INTEGER(iwp) ::  pr             !:
323       INTEGER(iwp) ::  wind_component !:
324       
325#if defined( __nopointer )
326       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var !:
327#else
328       REAL(wp), DIMENSION(:,:,:), POINTER ::  var
329#endif
330
331
332       IF ( .NOT. sloping_surface )  THEN
333!
334!--       Normal case: horizontal surface
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_wp * (    &
337                        ( var(k,j,i)   - ref_state(k)   ) / ref_state(k)   +   &
338                        ( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1)     &
339                                                                          )
340          ENDDO
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_wp * ( ( 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_wp * ( ( 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 END MODULE buoyancy_mod
Note: See TracBrowser for help on using the repository browser.