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

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

update of GPL copyright

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