source: palm/tags/release-3.8/SOURCE/pres.f90 @ 4121

Last change on this file since 4121 was 710, checked in by raasch, 13 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 22.3 KB
Line 
1 SUBROUTINE pres
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: pres.f90 710 2011-03-30 09:45:27Z schwenkel $
11!
12! 709 2011-03-30 09:31:40Z raasch
13! formatting adjustments
14!
15! 707 2011-03-29 11:39:40Z raasch
16! Calculation of weighted average of p is now handled in the same way
17! regardless of the number of ghost layers (advection scheme),
18! multigrid and sor method are using p_loc for iterative advancements of
19! pressure,
20! localsum calculation modified for proper OpenMP reduction,
21! bc_lr/ns replaced by bc_lr/ns_cyc
22!
23! 693 2011-03-08 09:..:..Z raasch
24! bugfix: weighting coefficient added to ibm branch
25!
26! 680 2011-02-04 23:16:06Z gryschka
27! bugfix: collective_wait
28!
29! 675 2011-01-19 10:56:55Z suehring
30! Removed bugfix while copying tend.
31!
32! 673 2011-01-18 16:19:48Z suehring
33! Weighting coefficients added for right computation of the pressure during
34! Runge-Kutta substeps.
35!
36! 667 2010-12-23 12:06:00Z suehring/gryschka
37! New allocation of tend when ws-scheme and multigrid is used. This is due to
38! reasons of perforance of the data_exchange. The same is done with p after
39! poismg is called.
40! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng when no
41! multigrid is used. Calls of exchange_horiz are modified.
42! bugfix: After pressure correction no volume flow correction in case of
43! non-cyclic boundary conditions
44! (has to be done only before pressure correction)
45! Call of SOR routine is referenced with ddzu_pres.
46!
47! 622 2010-12-10 08:08:13Z raasch
48! optional barriers included in order to speed up collective operations
49!
50! 151 2008-03-07 13:42:18Z raasch
51! Bugfix in volume flow control for non-cyclic boundary conditions
52!
53! 106 2007-08-16 14:30:26Z raasch
54! Volume flow conservation added for the remaining three outflow boundaries
55!
56! 85 2007-05-11 09:35:14Z raasch
57! Division through dt_3d replaced by multiplication of the inverse.
58! For performance optimisation, this is done in the loop calculating the
59! divergence instead of using a seperate loop.
60!
61! 75 2007-03-22 09:54:05Z raasch
62! Volume flow control for non-cyclic boundary conditions added (currently only
63! for the north boundary!!), 2nd+3rd argument removed from exchange horiz,
64! mean vertical velocity is removed in case of Neumann boundary conditions
65! both at the bottom and the top
66!
67! RCS Log replace by Id keyword, revision history cleaned up
68!
69! Revision 1.25  2006/04/26 13:26:12  raasch
70! OpenMP optimization (+localsum, threadsum)
71!
72! Revision 1.1  1997/07/24 11:24:44  raasch
73! Initial revision
74!
75!
76! Description:
77! ------------
78! Compute the divergence of the provisional velocity field. Solve the Poisson
79! equation for the perturbation pressure. Compute the final velocities using
80! this perturbation pressure. Compute the remaining divergence.
81!------------------------------------------------------------------------------!
82
83    USE arrays_3d
84    USE constants
85    USE control_parameters
86    USE cpulog
87    USE grid_variables
88    USE indices
89    USE interfaces
90    USE pegrid
91    USE poisfft_mod
92    USE poisfft_hybrid_mod
93    USE statistics
94
95    IMPLICIT NONE
96
97    INTEGER ::  i, j, k, sr
98
99    REAL    ::  ddt_3d, localsum, threadsum, d_weight_pres
100
101    REAL, DIMENSION(1:2) ::  volume_flow_l, volume_flow_offset
102    REAL, DIMENSION(1:nzt) ::  w_l, w_l_l
103
104
105    CALL cpu_log( log_point(8), 'pres', 'start' )
106
107
108    ddt_3d = 1.0 / dt_3d
109    d_weight_pres = 1.0 / weight_pres(intermediate_timestep_count)
110
111!
112!-- Multigrid method expects array d to have one ghost layer.
113!--
114    IF ( psolver == 'multigrid' )  THEN
115     
116       DEALLOCATE( d )
117       ALLOCATE( d(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ) 
118
119!
120!--    Since p is later used to hold the weighted average of the substeps, it
121!--    cannot be used in the iterative solver. Therefore, its initial value is
122!--    stored on p_loc, which is then iteratively advanced in every substep.
123       IF ( intermediate_timestep_count == 1 )  THEN
124          DO  i = nxl-1, nxr+1
125             DO  j = nys-1, nyn+1
126                DO  k = nzb, nzt+1
127                   p_loc(k,j,i) = p(k,j,i)
128                ENDDO
129             ENDDO
130          ENDDO
131       ENDIF
132       
133    ELSEIF ( psolver == 'sor'  .AND.  intermediate_timestep_count == 1 )  THEN
134
135!
136!--    Since p is later used to hold the weighted average of the substeps, it
137!--    cannot be used in the iterative solver. Therefore, its initial value is
138!--    stored on p_loc, which is then iteratively advanced in every substep.
139       p_loc = p
140
141    ENDIF
142
143!
144!-- Conserve the volume flow at the outflow in case of non-cyclic lateral
145!-- boundary conditions
146!-- WARNING: so far, this conservation does not work at the left/south
147!--          boundary if the topography at the inflow differs from that at the
148!--          outflow! For this case, volume_flow_area needs adjustment!
149!
150!-- Left/right
151    IF ( conserve_volume_flow  .AND.  ( outflow_l .OR. outflow_r ) )  THEN
152
153       volume_flow(1)   = 0.0
154       volume_flow_l(1) = 0.0
155
156       IF ( outflow_l )  THEN
157          i = 0
158       ELSEIF ( outflow_r )  THEN
159          i = nx+1
160       ENDIF
161
162       DO  j = nys, nyn
163!
164!--       Sum up the volume flow through the south/north boundary
165          DO  k = nzb_2d(j,i)+1, nzt
166             volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k)
167          ENDDO
168       ENDDO
169
170#if defined( __parallel )   
171       IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
172       CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, &
173                           MPI_SUM, comm1dy, ierr )   
174#else
175       volume_flow = volume_flow_l 
176#endif
177       volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) ) &
178                               / volume_flow_area(1)
179
180       DO  j = nysg, nyng
181          DO  k = nzb_2d(j,i)+1, nzt
182             u(k,j,i) = u(k,j,i) + volume_flow_offset(1)
183          ENDDO
184       ENDDO
185
186    ENDIF
187
188!
189!-- South/north
190    IF ( conserve_volume_flow  .AND.  ( outflow_n .OR. outflow_s ) )  THEN
191
192       volume_flow(2)   = 0.0
193       volume_flow_l(2) = 0.0
194
195       IF ( outflow_s )  THEN
196          j = 0
197       ELSEIF ( outflow_n )  THEN
198          j = ny+1
199       ENDIF
200
201       DO  i = nxl, nxr
202!
203!--       Sum up the volume flow through the south/north boundary
204          DO  k = nzb_2d(j,i)+1, nzt
205             volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k)
206          ENDDO
207       ENDDO
208
209#if defined( __parallel )   
210       IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
211       CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL, &
212                           MPI_SUM, comm1dx, ierr )   
213#else
214       volume_flow = volume_flow_l 
215#endif
216       volume_flow_offset(2) = ( volume_flow_initial(2) - volume_flow(2) )    &
217                               / volume_flow_area(2)
218
219       DO  i = nxlg, nxrg
220          DO  k = nzb_v_inner(j,i)+1, nzt
221             v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
222          ENDDO
223       ENDDO
224
225    ENDIF
226
227!
228!-- Remove mean vertical velocity
229    IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
230       IF ( simulated_time > 0.0 )  THEN ! otherwise nzb_w_inner not yet known
231          w_l = 0.0;  w_l_l = 0.0
232          DO  i = nxl, nxr
233             DO  j = nys, nyn
234                DO  k = nzb_w_inner(j,i)+1, nzt
235                   w_l_l(k) = w_l_l(k) + w(k,j,i)
236                ENDDO
237             ENDDO
238          ENDDO
239#if defined( __parallel )   
240          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
241          CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, &
242                              comm2d, ierr )
243#else
244          w_l = w_l_l 
245#endif
246          DO  k = 1, nzt
247             w_l(k) = w_l(k) / ngp_2dh_outer(k,0)
248          ENDDO
249          DO  i = nxlg, nxrg
250             DO  j = nysg, nyng
251                DO  k = nzb_w_inner(j,i)+1, nzt
252                   w(k,j,i) = w(k,j,i) - w_l(k)
253                ENDDO
254             ENDDO
255          ENDDO
256       ENDIF
257    ENDIF
258
259!
260!-- Compute the divergence of the provisional velocity field.
261    CALL cpu_log( log_point_s(1), 'divergence', 'start' )
262
263    IF ( psolver == 'multigrid' )  THEN
264       !$OMP PARALLEL DO SCHEDULE( STATIC )
265       DO  i = nxl-1, nxr+1
266          DO  j = nys-1, nyn+1
267             DO  k = nzb, nzt+1
268                d(k,j,i) = 0.0
269             ENDDO
270          ENDDO
271       ENDDO
272    ELSE
273       !$OMP PARALLEL DO SCHEDULE( STATIC )
274       DO  i = nxl, nxra
275          DO  j = nys, nyna
276             DO  k = nzb+1, nzta
277                d(k,j,i) = 0.0
278             ENDDO
279          ENDDO
280       ENDDO
281    ENDIF
282
283    localsum  = 0.0
284    threadsum = 0.0
285
286#if defined( __ibm )
287    !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
288    !$OMP DO SCHEDULE( STATIC )
289    DO  i = nxl, nxr
290       DO  j = nys, nyn
291          DO  k = nzb_s_inner(j,i)+1, nzt
292             d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
293                          ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
294                          ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d      &
295                                                                * d_weight_pres 
296          ENDDO
297!
298!--       Additional pressure boundary condition at the bottom boundary for
299!--       inhomogeneous Prandtl layer heat fluxes and temperatures, respectively
300!--       dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0.
301!--       This condition must not be applied at the start of a run, because then
302!--       flow_statistics has not yet been called and thus sums = 0.
303          IF ( ibc_p_b == 2  .AND.  sums(nzb+1,4) /= 0.0 )  THEN
304             k = nzb_s_inner(j,i)
305             d(k+1,j,i) = d(k+1,j,i) + (                                     &
306                                         ( usws(j,i+1) - usws(j,i) ) * ddx   &
307                                       + ( vsws(j+1,i) - vsws(j,i) ) * ddy   &
308                                       - g * ( pt(k+1,j,i) - sums(k+1,4) ) / &
309                                         sums(k+1,4)                         &
310                                       ) * ddzw(k+1) * ddt_3d * d_weight_pres 
311          ENDIF
312
313!
314!--       Compute possible PE-sum of divergences for flow_statistics
315          DO  k = nzb_s_inner(j,i)+1, nzt
316             threadsum = threadsum + ABS( d(k,j,i) )
317          ENDDO
318
319       ENDDO
320    ENDDO
321
322    localsum = localsum + threadsum * dt_3d * &
323                          weight_pres(intermediate_timestep_count)
324
325    !$OMP END PARALLEL
326#else
327    IF ( ibc_p_b == 2 .AND. sums(nzb+1,4) /= 0.0 )  THEN
328       !$OMP PARALLEL PRIVATE (i,j,k)
329       !$OMP DO SCHEDULE( STATIC )
330       DO  i = nxl, nxr
331          DO  j = nys, nyn
332             DO  k = nzb_s_inner(j,i)+1, nzt
333             d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
334                          ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
335                          ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d      &
336                                                                * d_weight_pres 
337             ENDDO
338          ENDDO
339!
340!--       Additional pressure boundary condition at the bottom boundary for
341!--       inhomogeneous Prandtl layer heat fluxes and temperatures, respectively
342!--       dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0.
343!--       This condition must not be applied at the start of a run, because then
344!--       flow_statistics has not yet been called and thus sums = 0.
345          DO  j = nys, nyn
346              k = nzb_s_inner(j,i)
347              d(k+1,j,i) = d(k+1,j,i) + (                        &
348                             ( usws(j,i+1) - usws(j,i) ) * ddx   &
349                           + ( vsws(j+1,i) - vsws(j,i) ) * ddy   &
350                           - g * ( pt(k+1,j,i) - sums(k+1,4) ) / &
351                             sums(k+1,4)                         &
352                                        ) * ddzw(k+1) * ddt_3d   &
353                                          * d_weight_pres 
354          ENDDO
355       ENDDO
356       !$OMP END PARALLEL
357
358    ELSE
359
360       !$OMP PARALLEL PRIVATE (i,j,k)
361       !$OMP DO SCHEDULE( STATIC )
362       DO  i = nxl, nxr
363          DO  j = nys, nyn
364             DO  k = nzb_s_inner(j,i)+1, nzt
365                d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
366                          ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
367                          ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d      &
368                                                                * d_weight_pres 
369             ENDDO
370          ENDDO
371       ENDDO
372       !$OMP END PARALLEL
373
374    ENDIF
375
376!
377!-- Compute possible PE-sum of divergences for flow_statistics
378    !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
379    !$OMP DO SCHEDULE( STATIC )
380    DO  i = nxl, nxr
381       DO  j = nys, nyn
382          DO  k = nzb+1, nzt
383             threadsum = threadsum + ABS( d(k,j,i) )
384          ENDDO
385       ENDDO
386    ENDDO
387    localsum = localsum + threadsum * dt_3d * &
388                          weight_pres(intermediate_timestep_count)
389    !$OMP END PARALLEL
390#endif
391
392!
393!-- For completeness, set the divergence sum of all statistic regions to those
394!-- of the total domain
395    sums_divold_l(0:statistic_regions) = localsum
396
397    CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
398
399!
400!-- Compute the pressure perturbation solving the Poisson equation
401    IF ( psolver(1:7) == 'poisfft' )  THEN
402
403!
404!--    Enlarge the size of tend, used as a working array for the transpositions
405       IF ( nxra > nxr  .OR.  nyna > nyn  .OR.  nza > nz )  THEN
406          DEALLOCATE( tend )
407          ALLOCATE( tend(1:nza,nys:nyna,nxl:nxra) )
408       ENDIF
409
410!
411!--    Solve Poisson equation via FFT and solution of tridiagonal matrices
412       IF ( psolver == 'poisfft' )  THEN
413!
414!--       Solver for 2d-decomposition
415          CALL poisfft( d, tend )
416       ELSEIF ( psolver == 'poisfft_hybrid' )  THEN 
417!
418!--       Solver for 1d-decomposition (using MPI and OpenMP).
419!--       The old hybrid-solver is still included here, as long as there
420!--       are some optimization problems in poisfft
421          CALL poisfft_hybrid( d )
422       ENDIF
423
424!
425!--    Resize tend to its normal size
426       IF ( nxra > nxr  .OR.  nyna > nyn  .OR.  nza > nz )  THEN
427          DEALLOCATE( tend )
428          ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
429       ENDIF
430
431!
432!--    Store computed perturbation pressure and set boundary condition in
433!--    z-direction
434       !$OMP PARALLEL DO
435       DO  i = nxl, nxr
436          DO  j = nys, nyn
437             DO  k = nzb+1, nzt
438                tend(k,j,i) = d(k,j,i)
439             ENDDO
440          ENDDO
441       ENDDO
442
443!
444!--    Bottom boundary:
445!--    This condition is only required for internal output. The pressure
446!--    gradient (dp(nzb+1)-dp(nzb))/dz is not used anywhere else.
447       IF ( ibc_p_b == 1 )  THEN
448!
449!--       Neumann (dp/dz = 0)
450          !$OMP PARALLEL DO
451          DO  i = nxlg, nxrg
452             DO  j = nysg, nyng
453                tend(nzb_s_inner(j,i),j,i) = tend(nzb_s_inner(j,i)+1,j,i)
454             ENDDO
455          ENDDO
456
457       ELSEIF ( ibc_p_b == 2 )  THEN
458!
459!--       Neumann condition for inhomogeneous surfaces,
460!--       here currently still in the form of a zero gradient. Actually
461!--       dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0 would have to be used for
462!--       the computation (cf. above: computation of divergences).
463          !$OMP PARALLEL DO
464          DO  i = nxlg, nxrg
465             DO  j = nysg, nyng
466                tend(nzb_s_inner(j,i),j,i) = tend(nzb_s_inner(j,i)+1,j,i)
467             ENDDO
468          ENDDO
469
470       ELSE
471!
472!--       Dirichlet
473          !$OMP PARALLEL DO
474          DO  i = nxlg, nxrg
475             DO  j = nysg, nyng
476                tend(nzb_s_inner(j,i),j,i) = 0.0
477             ENDDO
478          ENDDO
479
480       ENDIF
481
482!
483!--    Top boundary
484       IF ( ibc_p_t == 1 )  THEN
485!
486!--       Neumann
487          !$OMP PARALLEL DO
488          DO  i = nxlg, nxrg
489             DO  j = nysg, nyng
490                tend(nzt+1,j,i) = tend(nzt,j,i)
491             ENDDO
492          ENDDO
493
494       ELSE
495!
496!--       Dirichlet
497          !$OMP PARALLEL DO
498          DO  i = nxlg, nxrg
499             DO  j = nysg, nyng
500                tend(nzt+1,j,i) = 0.0
501             ENDDO
502          ENDDO
503
504       ENDIF
505
506!
507!--    Exchange boundaries for p
508       CALL exchange_horiz( tend, nbgp )
509     
510    ELSEIF ( psolver == 'sor' )  THEN
511
512!
513!--    Solve Poisson equation for perturbation pressure using SOR-Red/Black
514!--    scheme
515       CALL sor( d, ddzu_pres, ddzw, p_loc )
516       tend = p_loc
517
518    ELSEIF ( psolver == 'multigrid' )  THEN
519
520!
521!--    Solve Poisson equation for perturbation pressure using Multigrid scheme,
522!--    array tend is used to store the residuals, logical exchange_mg is used
523!--    to discern data exchange in multigrid ( 1 ghostpoint ) and normal grid
524!--    ( nbgp ghost points ).
525       CALL poismg( tend )
526
527!
528!--    Restore perturbation pressure on tend because this array is used
529!--    further below to correct the velocity fields
530       DO  i = nxl-1, nxr+1
531          DO  j = nys-1, nyn+1
532             DO  k = nzb, nzt+1
533                tend(k,j,i) = p_loc(k,j,i)
534             ENDDO
535          ENDDO
536       ENDDO
537
538    ENDIF
539
540!
541!-- Store perturbation pressure on array p, used for pressure data output.
542!-- Ghost layers are added in the output routines (except sor-method: see below)
543    IF ( intermediate_timestep_count == 1 )  THEN
544       !$OMP PARALLEL PRIVATE (i,j,k)
545       !$OMP DO
546       DO  i = nxl-1, nxr+1
547          DO  j = nys-1, nyn+1
548             DO  k = nzb, nzt+1
549                p(k,j,i) = tend(k,j,i) * &
550                           weight_substep(intermediate_timestep_count)
551             ENDDO
552          ENDDO
553       ENDDO
554       !$OMP END PARALLEL
555
556    ELSE 
557       !$OMP PARALLEL PRIVATE (i,j,k)
558       !$OMP DO
559       DO  i = nxl-1, nxr+1
560          DO  j = nys-1, nyn+1
561             DO  k = nzb, nzt+1
562                p(k,j,i) = p(k,j,i) + tend(k,j,i) * &
563                           weight_substep(intermediate_timestep_count)
564             ENDDO
565          ENDDO
566       ENDDO
567       !$OMP END PARALLEL
568
569    ENDIF
570       
571!
572!-- SOR-method needs ghost layers for the next timestep
573    IF ( psolver == 'sor' )  CALL exchange_horiz( p, nbgp )
574
575!
576!-- Correction of the provisional velocities with the current perturbation
577!-- pressure just computed
578    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc .OR. bc_ns_cyc ) )  THEN
579       volume_flow_l(1) = 0.0
580       volume_flow_l(2) = 0.0
581    ENDIF
582
583    !$OMP PARALLEL PRIVATE (i,j,k)
584    !$OMP DO
585    DO  i = nxl, nxr   
586       DO  j = nys, nyn
587          DO  k = nzb_w_inner(j,i)+1, nzt
588             w(k,j,i) = w(k,j,i) - dt_3d *                                 &
589                        ( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1) *      &
590                        weight_pres(intermediate_timestep_count)
591          ENDDO
592          DO  k = nzb_u_inner(j,i)+1, nzt
593             u(k,j,i) = u(k,j,i) - dt_3d *                                 &
594                        ( tend(k,j,i) - tend(k,j,i-1) ) * ddx *            &
595                        weight_pres(intermediate_timestep_count)
596          ENDDO
597          DO  k = nzb_v_inner(j,i)+1, nzt
598             v(k,j,i) = v(k,j,i) - dt_3d *                                 &
599                        ( tend(k,j,i) - tend(k,j-1,i) ) * ddy *            &
600                        weight_pres(intermediate_timestep_count)
601          ENDDO                                                         
602!
603!--       Sum up the volume flow through the right and north boundary
604          IF ( conserve_volume_flow  .AND.  bc_lr_cyc  .AND.  bc_ns_cyc  .AND. &
605               i == nx )  THEN
606             !$OMP CRITICAL
607             DO  k = nzb_2d(j,i) + 1, nzt
608                volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k)
609             ENDDO
610             !$OMP END CRITICAL
611          ENDIF
612          IF ( conserve_volume_flow  .AND.  bc_ns_cyc  .AND.  bc_lr_cyc  .AND. &
613               j == ny )  THEN
614             !$OMP CRITICAL
615             DO  k = nzb_2d(j,i) + 1, nzt
616                volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k)
617             ENDDO
618             !$OMP END CRITICAL
619          ENDIF
620
621       ENDDO
622    ENDDO
623    !$OMP END PARALLEL
624   
625!
626!-- Conserve the volume flow
627    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc  .AND.  bc_ns_cyc ) )  THEN
628
629#if defined( __parallel )   
630       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
631       CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 2, MPI_REAL, &
632                           MPI_SUM, comm2d, ierr ) 
633#else
634       volume_flow = volume_flow_l 
635#endif   
636
637       volume_flow_offset = ( volume_flow_initial - volume_flow ) / &
638                            volume_flow_area
639
640       !$OMP PARALLEL PRIVATE (i,j,k)
641       !$OMP DO
642       DO  i = nxl, nxr
643          DO  j = nys, nyn
644             DO  k = nzb_u_inner(j,i) + 1, nzt
645                u(k,j,i) = u(k,j,i) + volume_flow_offset(1)
646                v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
647             ENDDO
648          ENDDO
649       ENDDO
650
651       !$OMP END PARALLEL
652
653    ENDIF
654
655!
656!-- Exchange of boundaries for the velocities
657    CALL exchange_horiz( u, nbgp )
658    CALL exchange_horiz( v, nbgp )
659    CALL exchange_horiz( w, nbgp )
660
661!
662!-- Compute the divergence of the corrected velocity field,
663!-- a possible PE-sum is computed in flow_statistics
664    CALL cpu_log( log_point_s(1), 'divergence', 'start' )
665    sums_divnew_l = 0.0
666
667!
668!-- d must be reset to zero because it can contain nonzero values below the
669!-- topography
670    IF ( topography /= 'flat' )  d = 0.0
671
672    localsum  = 0.0
673    threadsum = 0.0
674
675    !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
676    !$OMP DO SCHEDULE( STATIC )
677#if defined( __ibm )
678    DO  i = nxl, nxr
679       DO  j = nys, nyn
680          DO  k = nzb_s_inner(j,i)+1, nzt
681             d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
682                        ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
683                        ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
684          ENDDO
685          DO  k = nzb+1, nzt
686             threadsum = threadsum + ABS( d(k,j,i) )
687          ENDDO
688       ENDDO
689    ENDDO
690#else
691    DO  i = nxl, nxr
692       DO  j = nys, nyn
693          DO  k = nzb_s_inner(j,i)+1, nzt
694             d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
695                        ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
696                        ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
697             threadsum = threadsum + ABS( d(k,j,i) )
698          ENDDO
699       ENDDO
700    ENDDO
701#endif
702
703    localsum = localsum + threadsum
704    !$OMP END PARALLEL
705
706!
707!-- For completeness, set the divergence sum of all statistic regions to those
708!-- of the total domain
709    sums_divnew_l(0:statistic_regions) = localsum
710
711    CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
712
713    CALL cpu_log( log_point(8), 'pres', 'stop' )
714   
715
716
717 END SUBROUTINE pres
Note: See TracBrowser for help on using the repository browser.