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

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

New:
---
Initial profiles can be used as reference state in the buoyancy term. New parameter
reference_state introduced. Calculation and handling of reference state in buoyancy term revised.
binary version for restart files changed from 3.9 to 3.9a (no downward compatibility!),
initial profile for rho added to hom (id=77)

Errors:


small bugfix for background communication (time_integration)

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