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

Last change on this file since 2956 was 2718, checked in by maronga, 7 years ago

deleting of deprecated files; headers updated where needed

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