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

Last change on this file since 1182 was 1182, checked in by raasch, 8 years ago

last commit documented, rc-file for example run updated

  • Property svn:keywords set to Id
File size: 14.0 KB
Line 
1 MODULE buoyancy_mod
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later 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-2012  Leibniz University Hannover
18!--------------------------------------------------------------------------------!
19!
20! Currrent revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: buoyancy.f90 1182 2013-06-14 09:07:24Z raasch $
27!
28! 1179 2013-06-14 05:57:58Z raasch
29! steering of reference state revised (var_reference and pr removed from
30! parameter list), use_reference renamed use_single_reference_value
31!
32! 1171 2013-05-30 11:27:45Z raasch
33! openacc statements added to use_reference-case in accelerator version
34!
35! 1153 2013-05-10 14:33:08Z raasch
36! code adjustments of accelerator version required by PGI 12.3 / CUDA 5.0
37!
38! 1128 2013-04-12 06:19:32Z raasch
39! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
40! j_north
41!
42! 1036 2012-10-22 13:43:42Z raasch
43! code put under GPL (PALM 3.9)
44!
45! 1015 2012-09-27 09:23:24Z raasch
46! accelerator version (*_acc) added
47!
48! 1010 2012-09-20 07:59:54Z raasch
49! cpp switch __nopointer added for pointer free version
50!
51! 622 2010-12-10 08:08:13Z raasch
52! optional barriers included in order to speed up collective operations
53!
54! 515 2010-03-18 02:30:38Z raasch
55! PGI-compiler creates SIGFPE in routine buoyancy, if opt>1 is used! Therefore,
56! opt=1 is forced by compiler-directive.
57!
58! 247 2009-02-27 14:01:30Z heinze
59! Output of messages replaced by message handling routine
60!
61! 132 2007-11-20 09:46:11Z letzel
62! Vertical scalar profiles now based on nzb_s_inner and ngp_2dh_s_inner.
63!
64! 106 2007-08-16 14:30:26Z raasch
65! i loop for u-component (sloping surface) is starting from nxlu (needed for
66! non-cyclic boundary conditions)
67!
68! 97 2007-06-21 08:23:15Z raasch
69! Routine reneralized to be used with temperature AND density:
70! argument theta renamed var, new argument var_reference,
71! use_pt_reference renamed use_reference,
72! calc_mean_pt_profile renamed calc_mean_profile
73!
74! 57 2007-03-09 12:05:41Z raasch
75! Reference temperature pt_reference can be used.
76!
77! RCS Log replace by Id keyword, revision history cleaned up
78!
79! Revision 1.19  2006/04/26 12:09:56  raasch
80! OpenMP optimization (one dimension added to sums_l)
81!
82! Revision 1.1  1997/08/29 08:56:48  raasch
83! Initial revision
84!
85!
86! Description:
87! ------------
88! Buoyancy term of the third component of the equation of motion.
89! WARNING: humidity is not regarded when using a sloping surface!
90!------------------------------------------------------------------------------!
91
92    PRIVATE
93    PUBLIC buoyancy, buoyancy_acc, calc_mean_profile
94
95    INTERFACE buoyancy
96       MODULE PROCEDURE buoyancy
97       MODULE PROCEDURE buoyancy_ij
98    END INTERFACE buoyancy
99
100    INTERFACE buoyancy_acc
101       MODULE PROCEDURE buoyancy_acc
102    END INTERFACE buoyancy_acc
103
104    INTERFACE calc_mean_profile
105       MODULE PROCEDURE calc_mean_profile
106    END INTERFACE calc_mean_profile
107
108 CONTAINS
109
110
111!------------------------------------------------------------------------------!
112! Call for all grid points
113!------------------------------------------------------------------------------!
114    SUBROUTINE buoyancy( var, wind_component )
115
116       USE arrays_3d
117       USE control_parameters
118       USE indices
119       USE pegrid
120
121       IMPLICIT NONE
122
123       INTEGER ::  i, j, k, wind_component
124#if defined( __nopointer )
125       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var
126#else
127       REAL, DIMENSION(:,:,:), POINTER ::  var
128#endif
129
130
131       IF ( .NOT. sloping_surface )  THEN
132!
133!--       Normal case: horizontal surface
134          DO  i = nxl, nxr
135             DO  j = nys, nyn
136                DO  k = nzb_s_inner(j,i)+1, nzt-1
137                   tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * (  &
138                          ( var(k,j,i)   - ref_state(k) )   / ref_state(k) +   &
139                          ( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1)   &
140                                                                            )
141                ENDDO
142             ENDDO
143          ENDDO
144
145       ELSE
146!
147!--       Buoyancy term for a surface with a slope in x-direction. The equations
148!--       for both the u and w velocity-component contain proportionate terms.
149!--       Temperature field at time t=0 serves as environmental temperature.
150!--       Reference temperature (pt_surface) is the one at the lower left corner
151!--       of the total domain.
152          IF ( wind_component == 1 )  THEN
153
154             DO  i = nxlu, nxr
155                DO  j = nys, nyn
156                   DO  k = nzb_s_inner(j,i)+1, nzt-1
157                      tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *      &
158                           0.5 * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
159                                 - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
160                                 ) / pt_surface
161                   ENDDO
162                ENDDO
163             ENDDO
164
165          ELSEIF ( wind_component == 3 )  THEN
166
167             DO  i = nxl, nxr
168                DO  j = nys, nyn
169                   DO  k = nzb_s_inner(j,i)+1, nzt-1
170                      tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *      &
171                           0.5 * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
172                                 - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
173                                 ) / pt_surface
174                   ENDDO
175                ENDDO
176            ENDDO
177
178          ELSE
179             
180             WRITE( message_string, * ) 'no term for component "',&
181                                       wind_component,'"'
182             CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
183
184          ENDIF
185
186       ENDIF
187
188    END SUBROUTINE buoyancy
189
190
191!------------------------------------------------------------------------------!
192! Call for all grid points - accelerator version
193!------------------------------------------------------------------------------!
194    SUBROUTINE buoyancy_acc( var, wind_component )
195
196       USE arrays_3d
197       USE control_parameters
198       USE indices
199       USE pegrid
200
201       IMPLICIT NONE
202
203       INTEGER ::  i, j, k, wind_component
204#if defined( __nopointer )
205       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var
206#else
207       REAL, DIMENSION(:,:,:), POINTER ::  var
208#endif
209
210
211       IF ( .NOT. sloping_surface )  THEN
212!
213!--       Normal case: horizontal surface
214          !$acc kernels present( nzb_s_inner, ref_state, tend, var )
215          !$acc loop
216          DO  i = i_left, i_right
217             DO  j = j_south, j_north
218                !$acc loop independent vector(32)
219                DO  k = nzb_s_inner(j,i)+1, nzt-1
220                   tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * (  &
221                          ( var(k,j,i)   - ref_state(k) )   / ref_state(k) +   &
222                          ( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1)   &
223                                                                            )
224                ENDDO
225             ENDDO
226          ENDDO
227          !$acc end kernels
228
229       ELSE
230!
231!--       Buoyancy term for a surface with a slope in x-direction. The equations
232!--       for both the u and w velocity-component contain proportionate terms.
233!--       Temperature field at time t=0 serves as environmental temperature.
234!--       Reference temperature (pt_surface) is the one at the lower left corner
235!--       of the total domain.
236          IF ( wind_component == 1 )  THEN
237
238             DO  i = nxlu, nxr
239                DO  j = nys, nyn
240                   DO  k = nzb_s_inner(j,i)+1, nzt-1
241                      tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *      &
242                           0.5 * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
243                                 - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
244                                 ) / pt_surface
245                   ENDDO
246                ENDDO
247             ENDDO
248
249          ELSEIF ( wind_component == 3 )  THEN
250
251             DO  i = nxl, nxr
252                DO  j = nys, nyn
253                   DO  k = nzb_s_inner(j,i)+1, nzt-1
254                      tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *      &
255                           0.5 * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
256                                 - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
257                                 ) / pt_surface
258                   ENDDO
259                ENDDO
260            ENDDO
261
262          ELSE
263
264             WRITE( message_string, * ) 'no term for component "',&
265                                       wind_component,'"'
266             CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
267
268          ENDIF
269
270       ENDIF
271
272    END SUBROUTINE buoyancy_acc
273
274
275!------------------------------------------------------------------------------!
276! Call for grid point i,j
277! ATTENTION: PGI-compiler creates SIGFPE if opt>1 is used! Therefore, opt=1 is
278!            forced by compiler-directive.
279!------------------------------------------------------------------------------!
280!pgi$r opt=1
281    SUBROUTINE buoyancy_ij( i, j, var, wind_component )
282
283       USE arrays_3d
284       USE control_parameters
285       USE indices
286       USE pegrid
287
288       IMPLICIT NONE
289
290       INTEGER ::  i, j, k, pr, wind_component
291#if defined( __nopointer )
292       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var
293#else
294       REAL, DIMENSION(:,:,:), POINTER ::  var
295#endif
296
297
298       IF ( .NOT. sloping_surface )  THEN
299!
300!--       Normal case: horizontal surface
301          DO  k = nzb_s_inner(j,i)+1, nzt-1
302              tend(k,j,i) = tend(k,j,i) + atmos_ocean_sign * g * 0.5 * (     &
303                        ( var(k,j,i)   - ref_state(k)   ) / ref_state(k)   + &
304                        ( var(k+1,j,i) - ref_state(k+1) ) / ref_state(k+1)   &
305                                                                       )
306          ENDDO
307
308       ELSE
309!
310!--       Buoyancy term for a surface with a slope in x-direction. The equations
311!--       for both the u and w velocity-component contain proportionate terms.
312!--       Temperature field at time t=0 serves as environmental temperature.
313!--       Reference temperature (pt_surface) is the one at the lower left corner
314!--       of the total domain.
315          IF ( wind_component == 1 )  THEN
316
317             DO  k = nzb_s_inner(j,i)+1, nzt-1
318                tend(k,j,i) = tend(k,j,i) + g * sin_alpha_surface *            &
319                           0.5 * ( ( pt(k,j,i-1)         + pt(k,j,i)         ) &
320                                 - ( pt_slope_ref(k,i-1) + pt_slope_ref(k,i) ) &
321                                 ) / pt_surface
322             ENDDO
323
324          ELSEIF ( wind_component == 3 )  THEN
325
326             DO  k = nzb_s_inner(j,i)+1, nzt-1
327                tend(k,j,i) = tend(k,j,i) + g * cos_alpha_surface *            &
328                           0.5 * ( ( pt(k,j,i)         + pt(k+1,j,i)         ) &
329                                 - ( pt_slope_ref(k,i) + pt_slope_ref(k+1,i) ) &
330                                 ) / pt_surface
331             ENDDO
332
333          ELSE
334
335             WRITE( message_string, * ) 'no term for component "',&
336                                       wind_component,'"'
337             CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
338
339          ENDIF
340
341       ENDIF
342
343    END SUBROUTINE buoyancy_ij
344
345
346    SUBROUTINE calc_mean_profile( var, pr )
347
348!------------------------------------------------------------------------------!
349! Description:
350! ------------
351! Calculate the horizontally averaged vertical temperature profile (pr=4 in case
352! of potential temperature, 44 in case of virtual potential temperature, and 64
353! in case of density (ocean runs)).
354!------------------------------------------------------------------------------!
355
356       USE arrays_3d,  ONLY: ref_state
357       USE control_parameters
358       USE indices
359       USE pegrid
360       USE statistics
361
362       IMPLICIT NONE
363
364       INTEGER ::  i, j, k, omp_get_thread_num, pr, tn
365#if defined( __nopointer )
366       REAL, DIMENSION(:,:,:) ::  var
367#else
368       REAL, DIMENSION(:,:,:), POINTER ::  var
369#endif
370
371!
372!--    Computation of the horizontally averaged profile of variable var, unless
373!--    already done by the relevant call from flow_statistics. The calculation
374!--    is done only for the first respective intermediate timestep in order to
375!--    spare communication time and to produce identical model results with jobs
376!--    which are calling flow_statistics at different time intervals.
377       IF ( .NOT. flow_statistics_called  .AND. &
378            intermediate_timestep_count == 1 )  THEN
379
380!
381!--       Horizontal average of variable var
382          tn           =   0  ! Default thread number in case of one thread
383          !$OMP PARALLEL PRIVATE( i, j, k, tn )
384!$        tn = omp_get_thread_num()
385          sums_l(:,pr,tn) = 0.0
386          !$OMP DO
387          DO  i = nxl, nxr
388             DO  j =  nys, nyn
389                DO  k = nzb_s_inner(j,i), nzt+1
390                   sums_l(k,pr,tn) = sums_l(k,pr,tn) + var(k,j,i)
391                ENDDO
392             ENDDO
393          ENDDO
394          !$OMP END PARALLEL
395
396          DO  i = 1, threads_per_task-1
397             sums_l(:,pr,0) = sums_l(:,pr,0) + sums_l(:,pr,i)
398          ENDDO
399
400#if defined( __parallel )
401
402          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
403          CALL MPI_ALLREDUCE( sums_l(nzb,pr,0), sums(nzb,pr), nzt+2-nzb, &
404                              MPI_REAL, MPI_SUM, comm2d, ierr )
405
406#else
407
408          sums(:,pr) = sums_l(:,pr,0)
409
410#endif
411
412          hom(:,1,pr,0) = sums(:,pr) / ngp_2dh_s_inner(:,0)
413
414       ENDIF
415
416       ref_state(:)  = hom(:,1,pr,0)   ! this is used in the buoyancy term
417
418    END SUBROUTINE calc_mean_profile
419
420 END MODULE buoyancy_mod
Note: See TracBrowser for help on using the repository browser.