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

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

last commit documented / copyright update

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