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

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

last commit documented

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