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

Last change on this file since 3274 was 3274, checked in by knoop, 3 years ago

Modularization of all bulk cloud physics code components

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