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

Last change on this file since 680 was 680, checked in by gryschka, 11 years ago

message string

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