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

Last change on this file since 3542 was 3538, checked in by suehring, 5 years ago

Remove unnecessary double masking of topography in advection and buoyancy terms

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