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

Last change on this file since 1365 was 1365, checked in by boeske, 7 years ago

large scale forcing enabled

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