source: palm/tags/release-3.2b/SOURCE/pres.f90 @ 388

Last change on this file since 388 was 90, checked in by raasch, 17 years ago

New:
---
Calculation and output of user-defined profiles. New &userpar parameters data_output_pr_user and max_pr_user.

check_parameters, flow_statistics, modules, parin, read_var_list, user_interface, write_var_list

Changed:


Division through dt_3d replaced by multiplication of the inverse. For performance optimisation, this is done in the loop calculating the divergence instead of using a seperate loop. (pres.f90) var_hom and var_sum renamed pr_palm.

data_output_profiles, flow_statistics, init_3d_model, modules, parin, pres, read_var_list, run_control, time_integration

Errors:


Bugfix: work_fft*_vec removed from some PRIVATE-declarations (poisfft).

Bugfix: field_chr renamed field_char (user_interface).

Bugfix: output of use_upstream_for_tke (header).

header, poisfft, user_interface

  • Property svn:keywords set to Id
File size: 17.3 KB
Line 
1 SUBROUTINE pres
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: pres.f90 90 2007-05-30 09:18:47Z raasch $
11!
12! 85 2007-05-11 09:35:14Z raasch
13! Division through dt_3d replaced by multiplication of the inverse.
14! For performance optimisation, this is done in the loop calculating the
15! divergence instead of using a seperate loop.
16!
17! 75 2007-03-22 09:54:05Z raasch
18! Volume flow control for non-cyclic boundary conditions added (currently only
19! for the north boundary!!), 2nd+3rd argument removed from exchange horiz,
20! mean vertical velocity is removed in case of Neumann boundary conditions
21! both at the bottom and the top
22!
23! RCS Log replace by Id keyword, revision history cleaned up
24!
25! Revision 1.25  2006/04/26 13:26:12  raasch
26! OpenMP optimization (+localsum, threadsum)
27!
28! Revision 1.1  1997/07/24 11:24:44  raasch
29! Initial revision
30!
31!
32! Description:
33! ------------
34! Compute the divergence of the provisional velocity field. Solve the Poisson
35! equation for the perturbation pressure. Compute the final velocities using
36! this perturbation pressure. Compute the remaining divergence.
37!------------------------------------------------------------------------------!
38
39    USE arrays_3d
40    USE constants
41    USE control_parameters
42    USE cpulog
43    USE grid_variables
44    USE indices
45    USE interfaces
46    USE pegrid
47    USE poisfft_mod
48    USE poisfft_hybrid_mod
49    USE statistics
50
51    IMPLICIT NONE
52
53    INTEGER ::  i, j, k, sr
54
55    REAL    ::  ddt_3d, localsum, threadsum
56
57    REAL, DIMENSION(1:2) ::  volume_flow_l, volume_flow_offset
58    REAL, DIMENSION(1:nzt) ::  w_l, w_l_l
59
60
61    CALL cpu_log( log_point(8), 'pres', 'start' )
62
63
64    ddt_3d = 1.0 / dt_3d
65
66!
67!-- Multigrid method needs additional grid points for the divergence array
68    IF ( psolver == 'multigrid' )  THEN
69       DEALLOCATE( d )
70       ALLOCATE( d(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
71    ENDIF
72
73!
74!-- Conserve the volume flow at the outflow in case of non-cyclic lateral
75!-- boundary conditions
76    IF ( conserve_volume_flow  .AND.  bc_ns == 'radiation/dirichlet')  THEN
77
78       volume_flow(2)   = 0.0
79       volume_flow_l(2) = 0.0
80
81       IF ( nyn == ny )  THEN
82          j = ny+1
83          DO  i = nxl, nxr
84!
85!--          Sum up the volume flow through the north boundary
86             DO  k = nzb_2d(j,i) + 1, nzt
87                volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzu(k)
88             ENDDO
89          ENDDO
90       ENDIF
91#if defined( __parallel )   
92       CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL, &
93                           MPI_SUM, comm1dx, ierr )   
94#else
95       volume_flow = volume_flow_l 
96#endif
97       volume_flow_offset(2) = ( volume_flow_initial(2) - volume_flow(2) )    &
98                               / volume_flow_area(2)                         
99
100       IF ( outflow_n )  THEN
101          j = nyn+1
102          DO  i = nxl, nxr
103             DO  k = nzb_v_inner(j,i) + 1, nzt
104                v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
105             ENDDO
106          ENDDO
107       ENDIF
108
109       CALL exchange_horiz( v )
110
111    ENDIF
112
113!
114!-- Remove mean vertical velocity
115    IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
116       IF ( simulated_time > 0.0 )  THEN ! otherwise nzb_w_inner is not yet known
117          w_l = 0.0;  w_l_l = 0.0
118          DO  i = nxl, nxr
119             DO  j = nys, nyn
120                DO  k = nzb_w_inner(j,i)+1, nzt
121                   w_l_l(k) = w_l_l(k) + w(k,j,i)
122                ENDDO
123             ENDDO
124          ENDDO
125#if defined( __parallel )   
126          CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, comm2d, &
127                              ierr )
128#else
129          w_l = w_l_l 
130#endif
131          DO  k = 1, nzt
132             w_l(k) = w_l(k) / ngp_2dh_outer(k,0)
133          ENDDO
134          DO  i = nxl-1, nxr+1
135             DO  j = nys-1, nyn+1
136                DO  k = nzb_w_inner(j,i)+1, nzt
137                   w(k,j,i) = w(k,j,i) - w_l(k)
138                ENDDO
139             ENDDO
140          ENDDO
141       ENDIF
142    ENDIF
143
144!
145!-- Compute the divergence of the provisional velocity field.
146    CALL cpu_log( log_point_s(1), 'divergence', 'start' )
147
148    IF ( psolver == 'multigrid' )  THEN
149       !$OMP PARALLEL DO SCHEDULE( STATIC )
150       DO  i = nxl-1, nxr+1
151          DO  j = nys-1, nyn+1
152             DO  k = nzb, nzt+1
153                d(k,j,i) = 0.0
154             ENDDO
155          ENDDO
156       ENDDO
157    ELSE
158       !$OMP PARALLEL DO SCHEDULE( STATIC )
159       DO  i = nxl, nxra
160          DO  j = nys, nyna
161             DO  k = nzb+1, nzta
162                d(k,j,i) = 0.0
163             ENDDO
164          ENDDO
165       ENDDO
166    ENDIF
167
168    localsum  = 0.0
169    threadsum = 0.0
170
171#if defined( __ibm )
172    !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
173    !$OMP DO SCHEDULE( STATIC )
174    DO  i = nxl, nxr
175       DO  j = nys, nyn
176          DO  k = nzb_s_inner(j,i)+1, nzt
177             d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
178                          ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
179                          ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d
180          ENDDO
181!
182!--       Additional pressure boundary condition at the bottom boundary for
183!--       inhomogeneous Prandtl layer heat fluxes and temperatures, respectively
184!--       dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0.
185!--       This condition must not be applied at the start of a run, because then
186!--       flow_statistics has not yet been called and thus sums = 0.
187          IF ( ibc_p_b == 2  .AND.  sums(nzb+1,4) /= 0.0 )  THEN
188             k = nzb_s_inner(j,i)
189             d(k+1,j,i) = d(k+1,j,i) + (                                     &
190                                         ( usws(j,i+1) - usws(j,i) ) * ddx   &
191                                       + ( vsws(j+1,i) - vsws(j,i) ) * ddy   &
192                                       - g * ( pt(k+1,j,i) - sums(k+1,4) ) / &
193                                         sums(k+1,4)                         &
194                                       ) * ddzw(k+1) * ddt_3d
195          ENDIF
196
197!
198!--       Compute possible PE-sum of divergences for flow_statistics
199          DO  k = nzb_s_inner(j,i)+1, nzt
200             threadsum = threadsum + ABS( d(k,j,i) )
201          ENDDO
202
203       ENDDO
204    ENDDO
205
206    localsum = ( localsum + threadsum ) * dt_3d
207    !$OMP END PARALLEL
208#else
209    IF ( ibc_p_b == 2 .AND. sums(nzb+1,4) /= 0.0 )  THEN
210       !$OMP PARALLEL PRIVATE (i,j,k)
211       !$OMP DO SCHEDULE( STATIC )
212       DO  i = nxl, nxr
213          DO  j = nys, nyn
214             DO  k = nzb_s_inner(j,i)+1, nzt
215                d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
216                             ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
217                             ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d
218             ENDDO
219          ENDDO
220!
221!--       Additional pressure boundary condition at the bottom boundary for
222!--       inhomogeneous Prandtl layer heat fluxes and temperatures, respectively
223!--       dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0.
224!--       This condition must not be applied at the start of a run, because then
225!--       flow_statistics has not yet been called and thus sums = 0.
226          DO  j = nys, nyn
227              k = nzb_s_inner(j,i)
228              d(k+1,j,i) = d(k+1,j,i) + (                        &
229                             ( usws(j,i+1) - usws(j,i) ) * ddx   &
230                           + ( vsws(j+1,i) - vsws(j,i) ) * ddy   &
231                           - g * ( pt(k+1,j,i) - sums(k+1,4) ) / &
232                             sums(k+1,4)                         &
233                                        ) * ddzw(k+1) * ddt_3d
234          ENDDO
235       ENDDO
236       !$OMP END PARALLEL
237
238    ELSE
239
240       !$OMP PARALLEL PRIVATE (i,j,k)
241       !$OMP DO SCHEDULE( STATIC )
242       DO  i = nxl, nxr
243          DO  j = nys, nyn
244             DO  k = nzb_s_inner(j,i)+1, nzt
245                d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
246                             ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
247                             ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d
248             ENDDO
249          ENDDO
250       ENDDO
251       !$OMP END PARALLEL
252
253    ENDIF
254
255!
256!-- Compute possible PE-sum of divergences for flow_statistics
257    !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
258    !$OMP DO SCHEDULE( STATIC )
259    DO  i = nxl, nxr
260       DO  j = nys, nyn
261          DO  k = nzb+1, nzt
262             threadsum = threadsum + ABS( d(k,j,i) )
263          ENDDO
264       ENDDO
265    ENDDO
266    localsum = ( localsum + threadsum ) * dt_3d
267    !$OMP END PARALLEL
268#endif
269
270!
271!-- For completeness, set the divergence sum of all statistic regions to those
272!-- of the total domain
273    sums_divold_l(0:statistic_regions) = localsum
274
275!
276!-- Determine absolute minimum/maximum (only for test cases, therefore as
277!-- comment line)
278!    CALL global_min_max( nzb+1, nzt, nys, nyn, nxl, nxr, d, 'abs', divmax, &
279!                         divmax_ijk )
280
281    CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
282
283!
284!-- Compute the pressure perturbation solving the Poisson equation
285    IF ( psolver(1:7) == 'poisfft' )  THEN
286
287!
288!--    Enlarge the size of tend, used as a working array for the transpositions
289       IF ( nxra > nxr  .OR.  nyna > nyn  .OR.  nza > nz )  THEN
290          DEALLOCATE( tend )
291          ALLOCATE( tend(1:nza,nys:nyna,nxl:nxra) )
292       ENDIF
293
294!
295!--    Solve Poisson equation via FFT and solution of tridiagonal matrices
296       IF ( psolver == 'poisfft' )  THEN
297!
298!--       Solver for 2d-decomposition
299          CALL poisfft( d, tend )
300       ELSEIF ( psolver == 'poisfft_hybrid' )  THEN 
301!
302!--       Solver for 1d-decomposition (using MPI and OpenMP).
303!--       The old hybrid-solver is still included here, as long as there
304!--       are some optimization problems in poisfft
305          CALL poisfft_hybrid( d )
306       ENDIF
307
308!
309!--    Resize tend to its normal size
310       IF ( nxra > nxr  .OR.  nyna > nyn  .OR.  nza > nz )  THEN
311          DEALLOCATE( tend )
312          ALLOCATE( tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
313       ENDIF
314
315!
316!--    Store computed perturbation pressure and set boundary condition in
317!--    z-direction
318       !$OMP PARALLEL DO
319       DO  i = nxl, nxr
320          DO  j = nys, nyn
321             DO  k = nzb+1, nzt
322                tend(k,j,i) = d(k,j,i)
323             ENDDO
324          ENDDO
325       ENDDO
326
327!
328!--    Bottom boundary:
329!--    This condition is only required for internal output. The pressure
330!--    gradient (dp(nzb+1)-dp(nzb))/dz is not used anywhere else.
331       IF ( ibc_p_b == 1 )  THEN
332!
333!--       Neumann (dp/dz = 0)
334          !$OMP PARALLEL DO
335          DO  i = nxl-1, nxr+1
336             DO  j = nys-1, nyn+1
337                tend(nzb_s_inner(j,i),j,i) = tend(nzb_s_inner(j,i)+1,j,i)
338             ENDDO
339          ENDDO
340
341       ELSEIF ( ibc_p_b == 2 )  THEN
342!
343!--       Neumann condition for inhomogeneous surfaces,
344!--       here currently still in the form of a zero gradient. Actually
345!--       dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0 would have to be used for
346!--       the computation (cf. above: computation of divergences).
347          !$OMP PARALLEL DO
348          DO  i = nxl-1, nxr+1
349             DO  j = nys-1, nyn+1
350                tend(nzb_s_inner(j,i),j,i) = tend(nzb_s_inner(j,i)+1,j,i)
351             ENDDO
352          ENDDO
353
354       ELSE
355!
356!--       Dirichlet
357          !$OMP PARALLEL DO
358          DO  i = nxl-1, nxr+1
359             DO  j = nys-1, nyn+1
360                tend(nzb_s_inner(j,i),j,i) = 0.0
361             ENDDO
362          ENDDO
363
364       ENDIF
365
366!
367!--    Top boundary
368       IF ( ibc_p_t == 1 )  THEN
369!
370!--       Neumann
371          !$OMP PARALLEL DO
372          DO  i = nxl-1, nxr+1
373             DO  j = nys-1, nyn+1
374                tend(nzt+1,j,i) = tend(nzt,j,i)
375             ENDDO
376          ENDDO
377
378       ELSE
379!
380!--       Dirichlet
381          !$OMP PARALLEL DO
382          DO  i = nxl-1, nxr+1
383             DO  j = nys-1, nyn+1
384                tend(nzt+1,j,i) = 0.0
385             ENDDO
386          ENDDO
387
388       ENDIF
389
390!
391!--    Exchange boundaries for p
392       CALL exchange_horiz( tend )
393     
394    ELSEIF ( psolver == 'sor' )  THEN
395
396!
397!--    Solve Poisson equation for perturbation pressure using SOR-Red/Black
398!--    scheme
399       CALL sor( d, ddzu, ddzw, p )
400       tend = p
401
402    ELSEIF ( psolver == 'multigrid' )  THEN
403
404!
405!--    Solve Poisson equation for perturbation pressure using Multigrid scheme,
406!--    array tend is used to store the residuals
407       CALL poismg( tend )
408 
409!
410!--    Restore perturbation pressure on tend because this array is used
411!--    further below to correct the velocity fields
412       tend = p
413
414    ENDIF
415
416!
417!-- Store perturbation pressure on array p, used in the momentum equations
418    IF ( psolver(1:7) == 'poisfft' )  THEN
419!
420!--    Here, only the values from the left and right boundaries are copied
421!--    The remaining values are copied in the following loop due to speed
422!--    optimization
423       !$OMP PARALLEL DO
424       DO  j = nys-1, nyn+1
425          DO  k = nzb, nzt+1
426             p(k,j,nxl-1) = tend(k,j,nxl-1)
427             p(k,j,nxr+1) = tend(k,j,nxr+1)
428          ENDDO
429       ENDDO
430    ENDIF
431
432!
433!-- Correction of the provisional velocities with the current perturbation
434!-- pressure just computed
435    IF ( conserve_volume_flow  .AND. &
436         ( bc_lr == 'cyclic'  .OR.  bc_ns == 'cyclic' ) )  THEN
437       volume_flow_l(1) = 0.0
438       volume_flow_l(2) = 0.0
439    ENDIF
440    !$OMP PARALLEL PRIVATE (i,j,k)
441    !$OMP DO
442    DO  i = nxl, nxr
443       IF ( psolver(1:7) == 'poisfft' )  THEN
444          DO  j = nys-1, nyn+1
445             DO  k = nzb, nzt+1
446                p(k,j,i) = tend(k,j,i)
447             ENDDO
448          ENDDO
449       ENDIF
450       DO  j = nys, nyn
451          DO  k = nzb_w_inner(j,i)+1, nzt
452             w(k,j,i) = w(k,j,i) - dt_3d * &
453                                   ( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1)
454          ENDDO
455          DO  k = nzb_u_inner(j,i)+1, nzt
456             u(k,j,i) = u(k,j,i) - dt_3d * ( tend(k,j,i) - tend(k,j,i-1) ) * ddx
457          ENDDO
458          DO  k = nzb_v_inner(j,i)+1, nzt
459             v(k,j,i) = v(k,j,i) - dt_3d * ( tend(k,j,i) - tend(k,j-1,i) ) * ddy
460          ENDDO
461
462!
463!--       Sum up the volume flow through the right and north boundary
464          IF ( conserve_volume_flow  .AND.  bc_lr == 'cyclic'  .AND. &
465               i == nx )  THEN
466             !$OMP CRITICAL
467             DO  k = nzb_2d(j,i) + 1, nzt
468                volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzu(k)
469             ENDDO
470             !$OMP END CRITICAL
471          ENDIF
472          IF ( conserve_volume_flow  .AND.  bc_ns == 'cyclic'  .AND. &
473               j == ny )  THEN
474             !$OMP CRITICAL
475             DO  k = nzb_2d(j,i) + 1, nzt
476                volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzu(k)
477             ENDDO
478             !$OMP END CRITICAL
479          ENDIF
480
481       ENDDO
482    ENDDO
483    !$OMP END PARALLEL
484
485!
486!-- Conserve the volume flow
487    IF ( conserve_volume_flow  .AND. &
488         ( bc_lr == 'cyclic'  .OR.  bc_ns == 'cyclic' ) )  THEN
489
490#if defined( __parallel )   
491       CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 2, MPI_REAL, &
492                           MPI_SUM, comm2d, ierr ) 
493#else
494       volume_flow = volume_flow_l 
495#endif   
496
497       volume_flow_offset = ( volume_flow_initial - volume_flow ) / &
498                            volume_flow_area
499
500       !$OMP PARALLEL PRIVATE (i,j,k)
501       !$OMP DO
502       DO  i = nxl, nxr
503          DO  j = nys, nyn
504             IF ( bc_lr == 'cyclic' )  THEN
505                DO  k = nzb_u_inner(j,i) + 1, nzt
506                   u(k,j,i) = u(k,j,i) + volume_flow_offset(1)
507                ENDDO
508             ENDIF
509             IF ( bc_ns == 'cyclic' )  THEN
510                DO  k = nzb_v_inner(j,i) + 1, nzt
511                   v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
512                ENDDO
513             ENDIF
514          ENDDO
515       ENDDO
516       !$OMP END PARALLEL
517
518    ENDIF
519
520!
521!-- Exchange of boundaries for the velocities
522    CALL exchange_horiz( u )
523    CALL exchange_horiz( v )
524    CALL exchange_horiz( w )
525
526!
527!-- Compute the divergence of the corrected velocity field,
528!-- a possible PE-sum is computed in flow_statistics
529    CALL cpu_log( log_point_s(1), 'divergence', 'start' )
530    sums_divnew_l = 0.0
531
532!
533!-- d must be reset to zero because it can contain nonzero values below the
534!-- topography
535    IF ( topography /= 'flat' )  d = 0.0
536
537    localsum  = 0.0
538    threadsum = 0.0
539
540    !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
541    !$OMP DO SCHEDULE( STATIC )
542#if defined( __ibm )
543    DO  i = nxl, nxr
544       DO  j = nys, nyn
545          DO  k = nzb_s_inner(j,i)+1, nzt
546             d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
547                        ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
548                        ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
549          ENDDO
550          DO  k = nzb+1, nzt
551             threadsum = threadsum + ABS( d(k,j,i) )
552          ENDDO
553       ENDDO
554    ENDDO
555#else
556    DO  i = nxl, nxr
557       DO  j = nys, nyn
558          DO  k = nzb_s_inner(j,i)+1, nzt
559             d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
560                        ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
561                        ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
562             threadsum = threadsum + ABS( d(k,j,i) )
563          ENDDO
564       ENDDO
565    ENDDO
566#endif
567    localsum = localsum + threadsum
568    !$OMP END PARALLEL
569
570!
571!-- For completeness, set the divergence sum of all statistic regions to those
572!-- of the total domain
573    sums_divnew_l(0:statistic_regions) = localsum
574
575    CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
576
577    CALL cpu_log( log_point(8), 'pres', 'stop' )
578
579
580 END SUBROUTINE pres
Note: See TracBrowser for help on using the repository browser.