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

Last change on this file since 1242 was 1242, checked in by heinze, 8 years ago

Last commmit documented

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