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

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

last commit documented

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