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

Last change on this file since 1132 was 1132, checked in by raasch, 9 years ago

r1028 documented

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