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

Last change on this file since 719 was 719, checked in by gryschka, 10 years ago

Bugfix in volume flow control for double cyclic boundary conditions

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