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

Last change on this file since 1172 was 1172, checked in by raasch, 11 years ago

last commit documented

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