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

Last change on this file since 778 was 778, checked in by fricke, 12 years ago

Modifications of the multigrid pressure solver

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