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

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

last commit documented

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