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

Last change on this file since 1179 was 1179, checked in by raasch, 8 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
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! steering of reference state revised (var_reference and pr removed from
23! parameter list), use_reference renamed use_single_reference_value
24!
25! Former revisions:
26! -----------------
27! $Id: buoyancy.f90 1179 2013-06-14 05:57:58Z raasch $
28!
29! 1171 2013-05-30 11:27:45Z raasch
30! openacc statements added to use_reference-case in accelerator version
31!
32! 1153 2013-05-10 14:33:08Z raasch
33! code adjustments of accelerator version required by PGI 12.3 / CUDA 5.0
34!
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!
39! 1036 2012-10-22 13:43:42Z raasch
40! code put under GPL (PALM 3.9)
41!
42! 1015 2012-09-27 09:23:24Z raasch
43! accelerator version (*_acc) added
44!
45! 1010 2012-09-20 07:59:54Z raasch
46! cpp switch __nopointer added for pointer free version
47!
48! 622 2010-12-10 08:08:13Z raasch
49! optional barriers included in order to speed up collective operations
50!
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!
55! 247 2009-02-27 14:01:30Z heinze
56! Output of messages replaced by message handling routine
57!
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!
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!
65! 97 2007-06-21 08:23:15Z raasch
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,
69! calc_mean_pt_profile renamed calc_mean_profile
70!
71! 57 2007-03-09 12:05:41Z raasch
72! Reference temperature pt_reference can be used.
73!
74! RCS Log replace by Id keyword, revision history cleaned up
75!
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
90    PUBLIC buoyancy, buoyancy_acc, calc_mean_profile
91
92    INTERFACE buoyancy
93       MODULE PROCEDURE buoyancy
94       MODULE PROCEDURE buoyancy_ij
95    END INTERFACE buoyancy
96
97    INTERFACE buoyancy_acc
98       MODULE PROCEDURE buoyancy_acc
99    END INTERFACE buoyancy_acc
100
101    INTERFACE calc_mean_profile
102       MODULE PROCEDURE calc_mean_profile
103    END INTERFACE calc_mean_profile
104
105 CONTAINS
106
107
108!------------------------------------------------------------------------------!
109! Call for all grid points
110!------------------------------------------------------------------------------!
111    SUBROUTINE buoyancy( var, wind_component )
112
113       USE arrays_3d
114       USE control_parameters
115       USE indices
116       USE pegrid
117
118       IMPLICIT NONE
119
120       INTEGER ::  i, j, k, wind_component
121#if defined( __nopointer )
122       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var
123#else
124       REAL, DIMENSION(:,:,:), POINTER ::  var
125#endif
126
127
128       IF ( .NOT. sloping_surface )  THEN
129!
130!--       Normal case: horizontal surface
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                                                                            )
138                ENDDO
139             ENDDO
140          ENDDO
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
151             DO  i = nxlu, 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 * 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
176             
177             WRITE( message_string, * ) 'no term for component "',&
178                                       wind_component,'"'
179             CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
180
181          ENDIF
182
183       ENDIF
184
185    END SUBROUTINE buoyancy
186
187
188!------------------------------------------------------------------------------!
189! Call for all grid points - accelerator version
190!------------------------------------------------------------------------------!
191    SUBROUTINE buoyancy_acc( var, wind_component )
192
193       USE arrays_3d
194       USE control_parameters
195       USE indices
196       USE pegrid
197
198       IMPLICIT NONE
199
200       INTEGER ::  i, j, k, wind_component
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
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                                                                            )
221                ENDDO
222             ENDDO
223          ENDDO
224          !$acc end kernels
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!------------------------------------------------------------------------------!
273! Call for grid point i,j
274! ATTENTION: PGI-compiler creates SIGFPE if opt>1 is used! Therefore, opt=1 is
275!            forced by compiler-directive.
276!------------------------------------------------------------------------------!
277!pgi$r opt=1
278    SUBROUTINE buoyancy_ij( i, j, var, wind_component )
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
288#if defined( __nopointer )
289       REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  var
290#else
291       REAL, DIMENSION(:,:,:), POINTER ::  var
292#endif
293
294
295       IF ( .NOT. sloping_surface )  THEN
296!
297!--       Normal case: horizontal surface
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
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
332             WRITE( message_string, * ) 'no term for component "',&
333                                       wind_component,'"'
334             CALL message( 'buoyancy', 'PA0159', 1, 2, 0, 6, 0 )
335
336          ENDIF
337
338       ENDIF
339
340    END SUBROUTINE buoyancy_ij
341
342
343    SUBROUTINE calc_mean_profile( var, pr )
344
345!------------------------------------------------------------------------------!
346! Description:
347! ------------
348! Calculate the horizontally averaged vertical temperature profile (pr=4 in case
349! of potential temperature, 44 in case of virtual potential temperature, and 64
350! in case of density (ocean runs)).
351!------------------------------------------------------------------------------!
352
353       USE arrays_3d,  ONLY: ref_state
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
362#if defined( __nopointer )
363       REAL, DIMENSION(:,:,:) ::  var
364#else
365       REAL, DIMENSION(:,:,:), POINTER ::  var
366#endif
367
368!
369!--    Computation of the horizontally averaged profile of variable var, unless
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!
378!--       Horizontal average of variable var
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
386                DO  k = nzb_s_inner(j,i), nzt+1
387                   sums_l(k,pr,tn) = sums_l(k,pr,tn) + var(k,j,i)
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
399          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
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
409          hom(:,1,pr,0) = sums(:,pr) / ngp_2dh_s_inner(:,0)
410
411       ENDIF
412
413       ref_state(:)  = hom(:,1,pr,0)   ! this is used in the buoyancy term
414
415    END SUBROUTINE calc_mean_profile
416
417 END MODULE buoyancy_mod
Note: See TracBrowser for help on using the repository browser.