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

Last change on this file since 684 was 682, checked in by gryschka, 13 years ago

removed print-statement for debugging

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