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

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

code has been put under the GNU General Public License (v3)

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