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

Last change on this file since 675 was 675, checked in by suehring, 11 years ago

Bugfix removed while copying p.

  • Property svn:keywords set to Id
File size: 23.2 KB
Line 
1 SUBROUTINE pres
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6! Removed bugfix while copying tend.
7!
8! Former revisions:
9! -----------------
10! $Id: pres.f90 675 2011-01-19 10:56:55Z suehring $
11!
12! 673 2011-01-18 16:19:48Z suehring
13! Weighting coefficients added for right computation of the pressure during
14! Runge-Kutta substeps.
15!
16! 667 2010-12-23 12:06:00Z suehring/gryschka
17! New allocation of tend when ws-scheme and multigrid is used. This is due to
18! reasons of perforance of the data_exchange. The same is done with p after
19! poismg is called.
20! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng when no
21! multigrid is used. Calls of exchange_horiz are modified.
22! bugfix: After pressure correction no volume flow correction in case of
23! non-cyclic boundary conditions
24! (has to be done only before pressure correction)
25! Call of SOR routine is referenced with ddzu_pres.
26!
27! 622 2010-12-10 08:08:13Z raasch
28! optional barriers included in order to speed up collective operations
29!
30! 151 2008-03-07 13:42:18Z raasch
31! Bugfix in volume flow control for non-cyclic boundary conditions
32!
33! 106 2007-08-16 14:30:26Z raasch
34! Volume flow conservation added for the remaining three outflow boundaries
35!
36! 85 2007-05-11 09:35:14Z raasch
37! Division through dt_3d replaced by multiplication of the inverse.
38! For performance optimisation, this is done in the loop calculating the
39! divergence instead of using a seperate loop.
40!
41! 75 2007-03-22 09:54:05Z raasch
42! Volume flow control for non-cyclic boundary conditions added (currently only
43! for the north boundary!!), 2nd+3rd argument removed from exchange horiz,
44! mean vertical velocity is removed in case of Neumann boundary conditions
45! both at the bottom and the top
46!
47! RCS Log replace by Id keyword, revision history cleaned up
48!
49! Revision 1.25  2006/04/26 13:26:12  raasch
50! OpenMP optimization (+localsum, threadsum)
51!
52! Revision 1.1  1997/07/24 11:24:44  raasch
53! Initial revision
54!
55!
56! Description:
57! ------------
58! Compute the divergence of the provisional velocity field. Solve the Poisson
59! equation for the perturbation pressure. Compute the final velocities using
60! this perturbation pressure. Compute the remaining divergence.
61!------------------------------------------------------------------------------!
62
63    USE arrays_3d
64    USE constants
65    USE control_parameters
66    USE cpulog
67    USE grid_variables
68    USE indices
69    USE interfaces
70    USE pegrid
71    USE poisfft_mod
72    USE poisfft_hybrid_mod
73    USE statistics
74
75    IMPLICIT NONE
76
77    INTEGER ::  i, j, k, sr
78
79    REAL    ::  ddt_3d, localsum, threadsum, d_weight_pres
80
81    REAL, DIMENSION(1:2) ::  volume_flow_l, volume_flow_offset
82    REAL, DIMENSION(1:nzt) ::  w_l, w_l_l
83
84
85    CALL cpu_log( log_point(8), 'pres', 'start' )
86
87
88    ddt_3d = 1.0 / dt_3d
89    d_weight_pres = 1. / weight_pres(intermediate_timestep_count)
90
91!
92!-- Multigrid method expects 1 additional grid point for the arrays
93!-- d, tend and p
94    IF ( psolver == 'multigrid' )  THEN
95     
96       DEALLOCATE( d )
97       ALLOCATE( d(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ) 
98       
99       IF ( ws_scheme_mom  .OR. ws_scheme_sca )  THEN
100       
101          DEALLOCATE( tend )
102          ALLOCATE( tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
103          DEALLOCATE( p )
104          ALLOCATE( p(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
105         
106       ENDIF
107       
108    ENDIF
109
110!
111!-- Conserve the volume flow at the outflow in case of non-cyclic lateral
112!-- boundary conditions
113!-- WARNING: so far, this conservation does not work at the left/south
114!--          boundary if the topography at the inflow differs from that at the
115!--          outflow! For this case, volume_flow_area needs adjustment!
116!
117!-- Left/right
118    IF ( conserve_volume_flow  .AND.  ( outflow_l  .OR. outflow_r ) )  THEN
119
120       volume_flow(1)   = 0.0
121       volume_flow_l(1) = 0.0
122
123       IF ( outflow_l )  THEN
124          i = 0
125       ELSEIF ( outflow_r )  THEN
126          i = nx+1
127       ENDIF
128
129       DO  j = nys, nyn
130!
131!--       Sum up the volume flow through the south/north boundary
132          DO  k = nzb_2d(j,i) + 1, nzt
133             volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k)
134          ENDDO
135       ENDDO
136
137#if defined( __parallel )   
138       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
139       CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, &
140                           MPI_SUM, comm1dy, ierr )   
141#else
142       volume_flow = volume_flow_l 
143#endif
144       volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) )    &
145                               / volume_flow_area(1)
146
147       DO  j = nysg, nyng
148          DO  k = nzb_2d(j,i) + 1, nzt
149             u(k,j,i) = u(k,j,i) + volume_flow_offset(1)
150          ENDDO
151       ENDDO
152
153    ENDIF
154
155!
156!-- South/north
157    IF ( conserve_volume_flow  .AND.  ( outflow_n  .OR. outflow_s ) )  THEN
158
159       volume_flow(2)   = 0.0
160       volume_flow_l(2) = 0.0
161
162       IF ( outflow_s )  THEN
163          j = 0
164       ELSEIF ( outflow_n )  THEN
165          j = ny+1
166       ENDIF
167
168       DO  i = nxl, nxr
169!
170!--       Sum up the volume flow through the south/north boundary
171          DO  k = nzb_2d(j,i) + 1, nzt
172             volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k)
173          ENDDO
174       ENDDO
175
176#if defined( __parallel )   
177       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
178       CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL, &
179                           MPI_SUM, comm1dx, ierr )   
180#else
181       volume_flow = volume_flow_l 
182#endif
183       volume_flow_offset(2) = ( volume_flow_initial(2) - volume_flow(2) )    &
184                               / volume_flow_area(2)
185
186       DO  i = nxlg, nxrg
187          DO  k = nzb_v_inner(j,i) + 1, nzt
188             v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
189          ENDDO
190       ENDDO
191
192    ENDIF
193
194!
195!-- Remove mean vertical velocity
196    IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
197       IF ( simulated_time > 0.0 )  THEN ! otherwise nzb_w_inner is not yet known
198          w_l = 0.0;  w_l_l = 0.0
199          DO  i = nxl, nxr
200             DO  j = nys, nyn
201                DO  k = nzb_w_inner(j,i)+1, nzt
202                   w_l_l(k) = w_l_l(k) + w(k,j,i)
203                ENDDO
204             ENDDO
205          ENDDO
206#if defined( __parallel )   
207          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
208          CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, comm2d, &
209                              ierr )
210#else
211          w_l = w_l_l 
212#endif
213          DO  k = 1, nzt
214             w_l(k) = w_l(k) / ngp_2dh_outer(k,0)
215          ENDDO
216          DO  i = nxlg, nxrg
217             DO  j = nysg, nyng
218                DO  k = nzb_w_inner(j,i)+1, nzt
219                   w(k,j,i) = w(k,j,i) - w_l(k)
220                ENDDO
221             ENDDO
222          ENDDO
223       ENDIF
224    ENDIF
225
226!
227!-- Compute the divergence of the provisional velocity field.
228    CALL cpu_log( log_point_s(1), 'divergence', 'start' )
229
230    IF ( psolver == 'multigrid' )  THEN
231       !$OMP PARALLEL DO SCHEDULE( STATIC )
232       DO  i = nxl-1, nxr+1
233          DO  j = nys-1, nyn+1
234             DO  k = nzb, nzt+1
235                d(k,j,i) = 0.0
236             ENDDO
237          ENDDO
238       ENDDO
239    ELSE
240       !$OMP PARALLEL DO SCHEDULE( STATIC )
241       DO  i = nxl, nxra
242          DO  j = nys, nyna
243             DO  k = nzb+1, nzta
244                d(k,j,i) = 0.0
245             ENDDO
246          ENDDO
247       ENDDO
248    ENDIF
249
250    localsum  = 0.0
251    threadsum = 0.0
252
253#if defined( __ibm )
254    !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
255    !$OMP DO SCHEDULE( STATIC )
256    DO  i = nxl, nxr
257       DO  j = nys, nyn
258          DO  k = nzb_s_inner(j,i)+1, nzt
259             d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
260                          ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
261                          ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d      &
262                                                                * d_weight_pres
263          ENDDO
264!
265!--       Additional pressure boundary condition at the bottom boundary for
266!--       inhomogeneous Prandtl layer heat fluxes and temperatures, respectively
267!--       dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0.
268!--       This condition must not be applied at the start of a run, because then
269!--       flow_statistics has not yet been called and thus sums = 0.
270          IF ( ibc_p_b == 2  .AND.  sums(nzb+1,4) /= 0.0 )  THEN
271             k = nzb_s_inner(j,i)
272             d(k+1,j,i) = d(k+1,j,i) + (                                     &
273                                         ( usws(j,i+1) - usws(j,i) ) * ddx   &
274                                       + ( vsws(j+1,i) - vsws(j,i) ) * ddy   &
275                                       - g * ( pt(k+1,j,i) - sums(k+1,4) ) / &
276                                         sums(k+1,4)                         &
277                                       ) * ddzw(k+1) * ddt_3d * d_weight_pres
278          ENDIF
279
280!
281!--       Compute possible PE-sum of divergences for flow_statistics
282          DO  k = nzb_s_inner(j,i)+1, nzt
283             threadsum = threadsum + ABS( d(k,j,i) )
284          ENDDO
285
286       ENDDO
287    ENDDO
288
289    localsum = ( localsum + threadsum ) * dt_3d
290    !$OMP END PARALLEL
291#else
292    IF ( ibc_p_b == 2 .AND. sums(nzb+1,4) /= 0.0 )  THEN
293       !$OMP PARALLEL PRIVATE (i,j,k)
294       !$OMP DO SCHEDULE( STATIC )
295       DO  i = nxl, nxr
296          DO  j = nys, nyn
297             DO  k = nzb_s_inner(j,i)+1, nzt
298             d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
299                          ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
300                          ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d      &
301                                                                * d_weight_pres
302             ENDDO
303          ENDDO
304!
305!--       Additional pressure boundary condition at the bottom boundary for
306!--       inhomogeneous Prandtl layer heat fluxes and temperatures, respectively
307!--       dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0.
308!--       This condition must not be applied at the start of a run, because then
309!--       flow_statistics has not yet been called and thus sums = 0.
310          DO  j = nys, nyn
311              k = nzb_s_inner(j,i)
312              d(k+1,j,i) = d(k+1,j,i) + (                        &
313                             ( usws(j,i+1) - usws(j,i) ) * ddx   &
314                           + ( vsws(j+1,i) - vsws(j,i) ) * ddy   &
315                           - g * ( pt(k+1,j,i) - sums(k+1,4) ) / &
316                             sums(k+1,4)                         &
317                                        ) * ddzw(k+1) * ddt_3d   &
318                                          * d_weight_pres
319          ENDDO
320       ENDDO
321       !$OMP END PARALLEL
322
323    ELSE
324
325       !$OMP PARALLEL PRIVATE (i,j,k)
326       !$OMP DO SCHEDULE( STATIC )
327       DO  i = nxl, nxr
328          DO  j = nys, nyn
329             DO  k = nzb_s_inner(j,i)+1, nzt
330                d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
331                          ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
332                          ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d      &
333                                                                * d_weight_pres
334             ENDDO
335          ENDDO
336       ENDDO
337       !$OMP END PARALLEL
338
339    ENDIF
340
341!
342!-- Compute possible PE-sum of divergences for flow_statistics
343    !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
344    !$OMP DO SCHEDULE( STATIC )
345    DO  i = nxl, nxr
346       DO  j = nys, nyn
347          DO  k = nzb+1, nzt
348             threadsum = threadsum + ABS( d(k,j,i) )
349          ENDDO
350       ENDDO
351    ENDDO
352    localsum = ( localsum + threadsum ) * dt_3d                               &
353               * weight_pres(intermediate_timestep_count)
354    !$OMP END PARALLEL
355#endif
356
357!
358!-- For completeness, set the divergence sum of all statistic regions to those
359!-- of the total domain
360    sums_divold_l(0:statistic_regions) = localsum
361
362!
363!-- Determine absolute minimum/maximum (only for test cases, therefore as
364!-- comment line)
365!    CALL global_min_max( nzb+1, nzt, nys, nyn, nxl, nxr, d, 'abs', divmax, &
366!                        divmax_ijk )
367
368    CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
369
370!
371!-- Compute the pressure perturbation solving the Poisson equation
372    IF ( psolver(1:7) == 'poisfft' )  THEN
373
374!
375!--    Enlarge the size of tend, used as a working array for the transpositions
376       IF ( nxra > nxr  .OR.  nyna > nyn  .OR.  nza > nz )  THEN
377          DEALLOCATE( tend )
378          ALLOCATE( tend(1:nza,nys:nyna,nxl:nxra) )
379       ENDIF
380
381!
382!--    Solve Poisson equation via FFT and solution of tridiagonal matrices
383       IF ( psolver == 'poisfft' )  THEN
384!
385!--       Solver for 2d-decomposition
386          CALL poisfft( d, tend )
387       ELSEIF ( psolver == 'poisfft_hybrid' )  THEN 
388!
389!--       Solver for 1d-decomposition (using MPI and OpenMP).
390!--       The old hybrid-solver is still included here, as long as there
391!--       are some optimization problems in poisfft
392          CALL poisfft_hybrid( d )
393       ENDIF
394
395!
396!--    Resize tend to its normal size
397       IF ( nxra > nxr  .OR.  nyna > nyn  .OR.  nza > nz )  THEN
398          DEALLOCATE( tend )
399          ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
400       ENDIF
401
402!
403!--    Store computed perturbation pressure and set boundary condition in
404!--    z-direction
405       !$OMP PARALLEL DO
406       DO  i = nxl, nxr
407          DO  j = nys, nyn
408             DO  k = nzb+1, nzt
409                tend(k,j,i) = d(k,j,i)
410             ENDDO
411          ENDDO
412       ENDDO
413
414!
415!--    Bottom boundary:
416!--    This condition is only required for internal output. The pressure
417!--    gradient (dp(nzb+1)-dp(nzb))/dz is not used anywhere else.
418       IF ( ibc_p_b == 1 )  THEN
419!
420!--       Neumann (dp/dz = 0)
421          !$OMP PARALLEL DO
422          DO  i = nxlg, nxrg
423             DO  j = nysg, nyng
424                tend(nzb_s_inner(j,i),j,i) = tend(nzb_s_inner(j,i)+1,j,i)
425             ENDDO
426          ENDDO
427
428       ELSEIF ( ibc_p_b == 2 )  THEN
429!
430!--       Neumann condition for inhomogeneous surfaces,
431!--       here currently still in the form of a zero gradient. Actually
432!--       dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0 would have to be used for
433!--       the computation (cf. above: computation of divergences).
434          !$OMP PARALLEL DO
435          DO  i = nxlg, nxrg
436             DO  j = nysg, nyng
437                tend(nzb_s_inner(j,i),j,i) = tend(nzb_s_inner(j,i)+1,j,i)
438             ENDDO
439          ENDDO
440
441       ELSE
442!
443!--       Dirichlet
444          !$OMP PARALLEL DO
445          DO  i = nxlg, nxrg
446             DO  j = nysg, nyng
447                tend(nzb_s_inner(j,i),j,i) = 0.0
448             ENDDO
449          ENDDO
450
451       ENDIF
452
453!
454!--    Top boundary
455       IF ( ibc_p_t == 1 )  THEN
456!
457!--       Neumann
458          !$OMP PARALLEL DO
459          DO  i = nxlg, nxrg
460             DO  j = nysg, nyng
461                tend(nzt+1,j,i) = tend(nzt,j,i)
462             ENDDO
463          ENDDO
464
465       ELSE
466!
467!--       Dirichlet
468          !$OMP PARALLEL DO
469          DO  i = nxlg, nxrg
470             DO  j = nysg, nyng
471                tend(nzt+1,j,i) = 0.0
472             ENDDO
473          ENDDO
474
475       ENDIF
476
477!
478!--    Exchange boundaries for p
479       CALL exchange_horiz( tend, nbgp )
480     
481    ELSEIF ( psolver == 'sor' )  THEN
482
483!
484!--    Solve Poisson equation for perturbation pressure using SOR-Red/Black
485!--    scheme
486       CALL sor( d, ddzu_pres, ddzw, p )
487       tend = p
488
489    ELSEIF ( psolver == 'multigrid' )  THEN
490
491!
492!--    Solve Poisson equation for perturbation pressure using Multigrid scheme,
493!--    array tend is used to store the residuals, logical exchange_mg is used
494!--    to discern data exchange in multigrid ( 1 ghostpoint ) and normal grid
495!--    ( nbgp ghost points ).
496       exchange_mg = .TRUE.
497       CALL poismg( tend )
498       exchange_mg = .FALSE.
499!
500!--    Restore perturbation pressure on tend because this array is used
501!--    further below to correct the velocity fields
502
503       tend = p
504       IF( ws_scheme_mom .OR. ws_scheme_sca )  THEN
505!       
506!--       Allocate p to its normal size and restore pressure.         
507          DEALLOCATE( p )
508          ALLOCATE( p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
509         
510       ENDIF
511
512    ENDIF
513
514!
515!-- Store perturbation pressure on array p, used in the momentum equations
516    IF ( psolver(1:7) == 'poisfft' )  THEN
517
518       IF ( intermediate_timestep_count == 1 )  THEN
519          !$OMP PARALLEL PRIVATE (i,j,k)
520          !$OMP DO
521          DO  i = nxlg, nxrg
522             DO  j = nysg, nyng
523                DO  k = nzb, nzt+1
524                   p(k,j,i) = tend(k,j,i)                                     &
525                            * weight_substep(intermediate_timestep_count)
526                ENDDO
527             ENDDO
528          ENDDO
529          !$OMP END PARALLEL
530         
531       ELSE 
532          !$OMP PARALLEL PRIVATE (i,j,k)
533          !$OMP DO
534          DO  i = nxlg, nxrg
535             DO  j = nysg, nyng
536                DO  k = nzb, nzt+1
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.