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

Last change on this file since 1682 was 1682, checked in by knoop, 6 years ago

Code annotations made doxygen readable

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