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

Last change on this file since 1003 was 1003, checked in by raasch, 9 years ago

subdomains must have identical size, i.e. grid_matching = "match" not allowed any more
parameter grid_matching removed
some obsolete variables removed

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