source: palm/trunk/SOURCE/pres.f90 @ 709

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

formatting adjustments

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