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

Last change on this file since 1320 was 1320, checked in by raasch, 7 years ago

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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