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

Last change on this file since 2089 was 2001, checked in by knoop, 8 years ago

last commit documented

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