source: palm/tags/release-5.0/SOURCE/buoyancy.f90 @ 3452

Last change on this file since 3452 was 2696, checked in by kanani, 6 years ago

Merge of branch palm4u into trunk

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