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

Last change on this file since 1952 was 1874, checked in by maronga, 8 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 13.8 KB
Line 
1!> @file buoyancy.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2016 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! ------------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: buoyancy.f90 1874 2016-04-18 14:51:10Z suehring $
26!
27! 1873 2016-04-18 14:50:06Z maronga
28! Module renamed (removed _mod)
29!
30!
31! 1850 2016-04-08 13:29:27Z maronga
32! Module renamed
33!
34!
35! 1682 2015-10-07 23:56:08Z knoop
36! Code annotations made doxygen readable
37!
38! 1374 2014-04-25 12:55:07Z raasch
39! missing variables added to ONLY list
40!
41! 1365 2014-04-22 15:03:56Z boeske
42! Calculation of reference state in subroutine calc_mean_profile moved to
43! subroutine time_integration,
44! subroutine calc_mean_profile moved to new file calc_mean_profile.f90
45!
46! 1353 2014-04-08 15:21:23Z heinze
47! REAL constants provided with KIND-attribute
48!
49! 1320 2014-03-20 08:40:49Z raasch
50! ONLY-attribute added to USE-statements,
51! kind-parameters added to all INTEGER and REAL declaration statements,
52! kinds are defined in new module kinds,
53! revision history before 2012 removed,
54! comment fields (!:) to be used for variable explanations added to
55! all variable declaration statements
56!
57! 1257 2013-11-08 15:18:40Z raasch
58! vector length (32) removed from openacc clause
59!
60! 1241 2013-10-30 11:36:58Z heinze
61! Generalize calc_mean_profile for wider use: use additional steering
62! character loc
63!
64! 1179 2013-06-14 05:57:58Z raasch
65! steering of reference state revised (var_reference and pr removed from
66! parameter list), use_reference renamed use_single_reference_value
67!
68! 1171 2013-05-30 11:27:45Z raasch
69! openacc statements added to use_reference-case in accelerator version
70!
71! 1153 2013-05-10 14:33:08Z raasch
72! code adjustments of accelerator version required by PGI 12.3 / CUDA 5.0
73!
74! 1128 2013-04-12 06:19:32Z raasch
75! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
76! j_north
77!
78! 1036 2012-10-22 13:43:42Z raasch
79! code put under GPL (PALM 3.9)
80!
81! 1015 2012-09-27 09:23:24Z raasch
82! accelerator version (*_acc) added
83!
84! 1010 2012-09-20 07:59:54Z raasch
85! cpp switch __nopointer added for pointer free version
86!
87! Revision 1.1  1997/08/29 08:56:48  raasch
88! Initial revision
89!
90!
91! Description:
92! ------------
93!> Buoyancy term of the third component of the equation of motion.
94!> @attention Humidity is not regarded when using a sloping surface!
95!------------------------------------------------------------------------------!
96 MODULE buoyancy_mod
97 
98
99    PRIVATE
100    PUBLIC buoyancy, buoyancy_acc
101
102    INTERFACE buoyancy
103       MODULE PROCEDURE buoyancy
104       MODULE PROCEDURE buoyancy_ij
105    END INTERFACE buoyancy
106
107    INTERFACE buoyancy_acc
108       MODULE PROCEDURE buoyancy_acc
109    END INTERFACE buoyancy_acc
110
111 CONTAINS
112
113
114!------------------------------------------------------------------------------!
115! Description:
116! ------------
117!> Call for all grid points
118!------------------------------------------------------------------------------!
119    SUBROUTINE buoyancy( var, wind_component )
120
121       USE arrays_3d,                                                          &
122           ONLY:  pt, pt_slope_ref, ref_state, tend
123
124       USE control_parameters,                                                 &
125           ONLY:  atmos_ocean_sign, cos_alpha_surface, g, message_string,      &
126                  pt_surface, sin_alpha_surface, sloping_surface
127
128       USE indices,                                                            &
129           ONLY:  nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nzb,       &
130                  nzb_s_inner, nzt
131
132       USE kinds
133
134       USE pegrid
135
136
137       IMPLICIT NONE
138
139       INTEGER(iwp) ::  i              !<
140       INTEGER(iwp) ::  j              !<
141       INTEGER(iwp) ::  k              !<
142       INTEGER(iwp) ::  wind_component !<
143       
144#if defined( __nopointer )
145       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var !<
146#else
147       REAL(wp), DIMENSION(:,:,:), POINTER ::  var
148#endif
149
150
151       IF ( .NOT. sloping_surface )  THEN
152!
153!--       Normal case: horizontal surface
154          DO  i = nxl, nxr
155             DO  j = nys, nyn
156                DO  k = nzb_s_inner(j,i)+1, nzt-1
157                   tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5_wp * &
158                          (                                                    &
159                             ( var(k,j,i)   - ref_state(k) )   / ref_state(k) + &
160                             ( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1) &
161                          )
162                ENDDO
163             ENDDO
164          ENDDO
165
166       ELSE
167!
168!--       Buoyancy term for a surface with a slope in x-direction. The equations
169!--       for both the u and w velocity-component contain proportionate terms.
170!--       Temperature field at time t=0 serves as environmental temperature.
171!--       Reference temperature (pt_surface) is the one at the lower left corner
172!--       of the total domain.
173          IF ( wind_component == 1 )  THEN
174
175             DO  i = nxlu, nxr
176                DO  j = nys, nyn
177                   DO  k = nzb_s_inner(j,i)+1, nzt-1
178                      tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *      &
179                           0.5_wp * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
180                                    - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
181                                    ) / pt_surface
182                   ENDDO
183                ENDDO
184             ENDDO
185
186          ELSEIF ( wind_component == 3 )  THEN
187
188             DO  i = nxl, nxr
189                DO  j = nys, nyn
190                   DO  k = nzb_s_inner(j,i)+1, nzt-1
191                      tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *      &
192                           0.5_wp * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
193                                    - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
194                                    ) / pt_surface
195                   ENDDO
196                ENDDO
197            ENDDO
198
199          ELSE
200             
201             WRITE( message_string, * ) 'no term for component "',             &
202                                       wind_component,'"'
203             CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
204
205          ENDIF
206
207       ENDIF
208
209    END SUBROUTINE buoyancy
210
211
212!------------------------------------------------------------------------------!
213! Description:
214! ------------
215!> Call for all grid points - accelerator version
216!------------------------------------------------------------------------------!
217    SUBROUTINE buoyancy_acc( var, wind_component )
218
219       USE arrays_3d,                                                          &
220           ONLY:  pt, pt_slope_ref, ref_state, tend
221
222       USE control_parameters,                                                 &
223           ONLY:  atmos_ocean_sign, cos_alpha_surface, g, message_string,      &
224                  pt_surface, sin_alpha_surface, sloping_surface
225
226       USE indices,                                                            &
227           ONLY:  i_left, i_right, j_north, j_south, nxl, nxlg, nxlu, nxr,     &
228                  nxrg, nyn, nyng, nys, nysg, nzb, nzb_s_inner, nzt
229
230       USE kinds
231
232       USE pegrid
233
234
235       IMPLICIT NONE
236
237       INTEGER(iwp) ::  i              !<
238       INTEGER(iwp) ::  j              !<
239       INTEGER(iwp) ::  k              !<
240       INTEGER(iwp) ::  wind_component !<
241       
242#if defined( __nopointer )
243       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var !<
244#else
245       REAL(wp), DIMENSION(:,:,:), POINTER ::  var
246#endif
247
248
249       IF ( .NOT. sloping_surface )  THEN
250!
251!--       Normal case: horizontal surface
252          !$acc kernels present( nzb_s_inner, ref_state, tend, var )
253          !$acc loop
254          DO  i = i_left, i_right
255             DO  j = j_south, j_north
256                !$acc loop independent vector
257                DO  k = nzb_s_inner(j,i)+1, nzt-1
258                   tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5_wp * &
259                          (                                                    &
260                             ( var(k,j,i)   - ref_state(k) )   / ref_state(k) + &
261                             ( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1) &
262                          )
263                ENDDO
264             ENDDO
265          ENDDO
266          !$acc end kernels
267
268       ELSE
269!
270!--       Buoyancy term for a surface with a slope in x-direction. The equations
271!--       for both the u and w velocity-component contain proportionate terms.
272!--       Temperature field at time t=0 serves as environmental temperature.
273!--       Reference temperature (pt_surface) is the one at the lower left corner
274!--       of the total domain.
275          IF ( wind_component == 1 )  THEN
276
277             DO  i = nxlu, nxr
278                DO  j = nys, nyn
279                   DO  k = nzb_s_inner(j,i)+1, nzt-1
280                      tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *      &
281                           0.5_wp * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
282                                    - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
283                                    ) / pt_surface
284                   ENDDO
285                ENDDO
286             ENDDO
287
288          ELSEIF ( wind_component == 3 )  THEN
289
290             DO  i = nxl, nxr
291                DO  j = nys, nyn
292                   DO  k = nzb_s_inner(j,i)+1, nzt-1
293                      tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *      &
294                           0.5_wp * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
295                                    - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
296                                    ) / pt_surface
297                   ENDDO
298                ENDDO
299            ENDDO
300
301          ELSE
302
303             WRITE( message_string, * ) 'no term for component "',             &
304                                       wind_component,'"'
305             CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
306
307          ENDIF
308
309       ENDIF
310
311    END SUBROUTINE buoyancy_acc
312
313
314!------------------------------------------------------------------------------!
315! Description:
316! ------------
317!> Call for grid point i,j
318!> @attention PGI-compiler creates SIGFPE if opt>1 is used! Therefore, opt=1 is
319!>            forced by compiler-directive.
320!------------------------------------------------------------------------------!
321!pgi$r opt=1
322    SUBROUTINE buoyancy_ij( i, j, var, wind_component )
323
324       USE arrays_3d,                                                          &
325           ONLY:  pt, pt_slope_ref, ref_state, tend
326
327       USE control_parameters,                                                 &
328           ONLY:  atmos_ocean_sign, cos_alpha_surface, g, message_string,      &
329                  pt_surface, sin_alpha_surface, sloping_surface
330
331       USE indices,                                                            &
332           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzb_s_inner, nzt
333
334       USE kinds
335
336       USE pegrid
337
338
339       IMPLICIT NONE
340
341       INTEGER(iwp) ::  i              !<
342       INTEGER(iwp) ::  j              !<
343       INTEGER(iwp) ::  k              !<
344       INTEGER(iwp) ::  pr             !<
345       INTEGER(iwp) ::  wind_component !<
346       
347#if defined( __nopointer )
348       REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var !<
349#else
350       REAL(wp), DIMENSION(:,:,:), POINTER ::  var
351#endif
352
353
354       IF ( .NOT. sloping_surface )  THEN
355!
356!--       Normal case: horizontal surface
357          DO  k = nzb_s_inner(j,i)+1, nzt-1
358              tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5_wp * (    &
359                        ( var(k,j,i)   - ref_state(k)   ) / ref_state(k)   +   &
360                        ( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1)     &
361                                                                          )
362          ENDDO
363
364       ELSE
365!
366!--       Buoyancy term for a surface with a slope in x-direction. The equations
367!--       for both the u and w velocity-component contain proportionate terms.
368!--       Temperature field at time t=0 serves as environmental temperature.
369!--       Reference temperature (pt_surface) is the one at the lower left corner
370!--       of the total domain.
371          IF ( wind_component == 1 )  THEN
372
373             DO  k = nzb_s_inner(j,i)+1, nzt-1
374                tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *            &
375                           0.5_wp * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
376                                    - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
377                                    ) / pt_surface
378             ENDDO
379
380          ELSEIF ( wind_component == 3 )  THEN
381
382             DO  k = nzb_s_inner(j,i)+1, nzt-1
383                tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *            &
384                           0.5_wp * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
385                                    - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
386                                    ) / pt_surface
387             ENDDO
388
389          ELSE
390
391             WRITE( message_string, * ) 'no term for component "',             &
392                                       wind_component,'"'
393             CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
394
395          ENDIF
396
397       ENDIF
398
399    END SUBROUTINE buoyancy_ij
400
401 END MODULE buoyancy_mod
Note: See TracBrowser for help on using the repository browser.