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

Last change on this file since 1032 was 1017, checked in by raasch, 12 years ago

last commit documented

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