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

Last change on this file since 2000 was 2000, checked in by knoop, 5 years ago

Forced header and separation lines into 80 columns

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