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

Last change on this file since 693 was 693, checked in by raasch, 11 years ago

New:
---
adjustments for ibmkisti (pres, batch_scp, mbuild, mrun, subjob)

new default configuration files .mrun.config.yonsei2011 and .config_block_for_ibmkisti

Changed:


Errors:


bugfix: weighting coefficient added to ibm branch (pres)

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