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

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

last commit documented

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