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

Last change on this file since 708 was 708, checked in by raasch, 11 years ago

last commit documented

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