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

Last change on this file since 3661 was 3655, checked in by knoop, 5 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

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