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

Last change on this file since 1873 was 1873, checked in by maronga, 5 years ago

revised renaming of modules

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