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

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

last commit documented

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