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

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