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

Last change on this file since 1126 was 1037, checked in by raasch, 12 years ago

last commit documented

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