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

Last change on this file since 707 was 707, checked in by raasch, 10 years ago

New:
---

In case of multigrid method, on coarse grid levels, gathered data are
identically processed on all PEs (before, on PE0 only), so that the subsequent
scattering of data is not neccessary any more. (modules, init_pegrid, poismg)

Changed:


Calculation of weighted average of p is now handled in the same way
regardless of the number of ghost layers (advection scheme). (pres)

multigrid and sor method are using p_loc for iterative
advancements of pressure. p_sub removed. (init_3d_model, modules, poismg, pres, sor)

bc_lr and bc_ns replaced by bc_lr_dirrad, bc_lr_raddir, bc_ns_dirrad, bc_ns_raddir
for speed optimization. (calc_spectra, check_parameters, exchange_horiz,
exchange_horiz_2d, header, init_3d_model, init_grid, init_pegrid, modules,
poismg, pres, sor, time_integration, timestep)

grid_level directly used as index for MPI data type arrays. (exchange_horiz,
poismg)

initial assignments of zero to array p for iterative solvers only (init_3d_model)

Errors:


localsum calculation modified for proper OpenMP reduction. (pres)

Bugfix: bottom (nzb) and top (nzt+1) boundary conditions set in routines
resid and restrict. They were missed before, which may have led to
unpredictable results. (poismg)

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