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

Last change on this file since 2232 was 2232, checked in by suehring, 4 years ago

Adjustments according new topography and surface-modelling concept implemented

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