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

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

last commit documented

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