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

Last change on this file since 673 was 673, checked in by suehring, 13 years ago

Right computation of the pressure using Runge-Kutta weighting coefficients. Consideration of the pressure gradient during the time integration removed. Removed bugfix concerning velocity variances.

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