source: palm/trunk/SOURCE/buoyancy_mod.f90 @ 1850

Last change on this file since 1850 was 1850, checked in by maronga, 6 years ago

added _mod string to several filenames to meet the naming convection for modules

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