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

Last change on this file since 1011 was 1004, checked in by raasch, 12 years ago

last commit documented

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