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

Last change on this file since 695 was 695, checked in by raasch, 14 years ago

comment line changed

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