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

Last change on this file since 1414 was 1375, checked in by raasch, 11 years ago

last commit documented

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