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

Last change on this file since 779 was 779, checked in by fricke, 10 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 23.3 KB
Line 
1 SUBROUTINE pres
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: pres.f90 779 2011-11-07 14:30:01Z fricke $
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, nxra
282          DO  j = nys, nyna
283             DO  k = nzb+1, nzta
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!--    Enlarge the size of tend, used as a working array for the transpositions
412       IF ( nxra > nxr  .OR.  nyna > nyn  .OR.  nza > nz )  THEN
413          DEALLOCATE( tend )
414          ALLOCATE( tend(1:nza,nys:nyna,nxl:nxra) )
415       ENDIF
416
417!
418!--    Solve Poisson equation via FFT and solution of tridiagonal matrices
419       IF ( psolver == 'poisfft' )  THEN
420!
421!--       Solver for 2d-decomposition
422          CALL poisfft( d, tend )
423       ELSEIF ( psolver == 'poisfft_hybrid' )  THEN 
424!
425!--       Solver for 1d-decomposition (using MPI and OpenMP).
426!--       The old hybrid-solver is still included here, as long as there
427!--       are some optimization problems in poisfft
428          CALL poisfft_hybrid( d )
429       ENDIF
430
431!
432!--    Resize tend to its normal size
433       IF ( nxra > nxr  .OR.  nyna > nyn  .OR.  nza > nz )  THEN
434          DEALLOCATE( tend )
435          ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
436       ENDIF
437
438!
439!--    Store computed perturbation pressure and set boundary condition in
440!--    z-direction
441       !$OMP PARALLEL DO
442       DO  i = nxl, nxr
443          DO  j = nys, nyn
444             DO  k = nzb+1, nzt
445                tend(k,j,i) = d(k,j,i)
446             ENDDO
447          ENDDO
448       ENDDO
449
450!
451!--    Bottom boundary:
452!--    This condition is only required for internal output. The pressure
453!--    gradient (dp(nzb+1)-dp(nzb))/dz is not used anywhere else.
454       IF ( ibc_p_b == 1 )  THEN
455!
456!--       Neumann (dp/dz = 0)
457          !$OMP PARALLEL DO
458          DO  i = nxlg, nxrg
459             DO  j = nysg, nyng
460                tend(nzb_s_inner(j,i),j,i) = tend(nzb_s_inner(j,i)+1,j,i)
461             ENDDO
462          ENDDO
463
464       ELSEIF ( ibc_p_b == 2 )  THEN
465!
466!--       Neumann condition for inhomogeneous surfaces,
467!--       here currently still in the form of a zero gradient. Actually
468!--       dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0 would have to be used for
469!--       the computation (cf. above: computation of divergences).
470          !$OMP PARALLEL DO
471          DO  i = nxlg, nxrg
472             DO  j = nysg, nyng
473                tend(nzb_s_inner(j,i),j,i) = tend(nzb_s_inner(j,i)+1,j,i)
474             ENDDO
475          ENDDO
476
477       ELSE
478!
479!--       Dirichlet
480          !$OMP PARALLEL DO
481          DO  i = nxlg, nxrg
482             DO  j = nysg, nyng
483                tend(nzb_s_inner(j,i),j,i) = 0.0
484             ENDDO
485          ENDDO
486
487       ENDIF
488
489!
490!--    Top boundary
491       IF ( ibc_p_t == 1 )  THEN
492!
493!--       Neumann
494          !$OMP PARALLEL DO
495          DO  i = nxlg, nxrg
496             DO  j = nysg, nyng
497                tend(nzt+1,j,i) = tend(nzt,j,i)
498             ENDDO
499          ENDDO
500
501       ELSE
502!
503!--       Dirichlet
504          !$OMP PARALLEL DO
505          DO  i = nxlg, nxrg
506             DO  j = nysg, nyng
507                tend(nzt+1,j,i) = 0.0
508             ENDDO
509          ENDDO
510
511       ENDIF
512
513!
514!--    Exchange boundaries for p
515       CALL exchange_horiz( tend, nbgp )
516     
517    ELSEIF ( psolver == 'sor' )  THEN
518
519!
520!--    Solve Poisson equation for perturbation pressure using SOR-Red/Black
521!--    scheme
522       CALL sor( d, ddzu_pres, ddzw, p_loc )
523       tend = p_loc
524
525    ELSEIF ( psolver == 'multigrid' )  THEN
526
527!
528!--    Solve Poisson equation for perturbation pressure using Multigrid scheme,
529!--    array tend is used to store the residuals, logical exchange_mg is used
530!--    to discern data exchange in multigrid ( 1 ghostpoint ) and normal grid
531!--    ( nbgp ghost points ).
532
533!--    If the number of grid points of the gathered grid, which is collected
534!--    on PE0, is larger than the number of grid points of an PE, than array
535!--    tend will be enlarged.
536       IF ( gathered_size > subdomain_size )  THEN
537          DEALLOCATE( tend )
538          ALLOCATE( tend(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg(          &
539                    mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1,&
540                    nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg(                   &
541                    mg_switch_to_pe0_level)+1) )
542       ENDIF
543
544       CALL poismg( tend )
545
546       IF ( gathered_size > subdomain_size )  THEN
547          DEALLOCATE( tend )
548          ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
549       ENDIF
550
551!
552!--    Restore perturbation pressure on tend because this array is used
553!--    further below to correct the velocity fields
554       DO  i = nxl-1, nxr+1
555          DO  j = nys-1, nyn+1
556             DO  k = nzb, nzt+1
557                tend(k,j,i) = p_loc(k,j,i)
558             ENDDO
559          ENDDO
560       ENDDO
561
562    ENDIF
563
564!
565!-- Store perturbation pressure on array p, used for pressure data output.
566!-- Ghost layers are added in the output routines (except sor-method: see below)
567    IF ( intermediate_timestep_count == 1 )  THEN
568       !$OMP PARALLEL PRIVATE (i,j,k)
569       !$OMP DO
570       DO  i = nxl-1, nxr+1
571          DO  j = nys-1, nyn+1
572             DO  k = nzb, nzt+1
573                p(k,j,i) = tend(k,j,i) * &
574                           weight_substep(intermediate_timestep_count)
575             ENDDO
576          ENDDO
577       ENDDO
578       !$OMP END PARALLEL
579
580    ELSE 
581       !$OMP PARALLEL PRIVATE (i,j,k)
582       !$OMP DO
583       DO  i = nxl-1, nxr+1
584          DO  j = nys-1, nyn+1
585             DO  k = nzb, nzt+1
586                p(k,j,i) = p(k,j,i) + tend(k,j,i) * &
587                           weight_substep(intermediate_timestep_count)
588             ENDDO
589          ENDDO
590       ENDDO
591       !$OMP END PARALLEL
592
593    ENDIF
594       
595!
596!-- SOR-method needs ghost layers for the next timestep
597    IF ( psolver == 'sor' )  CALL exchange_horiz( p, nbgp )
598
599!
600!-- Correction of the provisional velocities with the current perturbation
601!-- pressure just computed
602    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc .OR. bc_ns_cyc ) )  THEN
603       volume_flow_l(1) = 0.0
604       volume_flow_l(2) = 0.0
605    ENDIF
606
607    !$OMP PARALLEL PRIVATE (i,j,k)
608    !$OMP DO
609    DO  i = nxl, nxr   
610       DO  j = nys, nyn
611          DO  k = nzb_w_inner(j,i)+1, nzt
612             w(k,j,i) = w(k,j,i) - dt_3d *                                 &
613                        ( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1) *      &
614                        weight_pres(intermediate_timestep_count)
615          ENDDO
616          DO  k = nzb_u_inner(j,i)+1, nzt
617             u(k,j,i) = u(k,j,i) - dt_3d *                                 &
618                        ( tend(k,j,i) - tend(k,j,i-1) ) * ddx *            &
619                        weight_pres(intermediate_timestep_count)
620          ENDDO
621          DO  k = nzb_v_inner(j,i)+1, nzt
622             v(k,j,i) = v(k,j,i) - dt_3d *                                 &
623                        ( tend(k,j,i) - tend(k,j-1,i) ) * ddy *            &
624                        weight_pres(intermediate_timestep_count)
625          ENDDO                                                         
626!
627!--       Sum up the volume flow through the right and north boundary
628          IF ( conserve_volume_flow  .AND.  bc_lr_cyc  .AND.  bc_ns_cyc  .AND. &
629               i == nx )  THEN
630             !$OMP CRITICAL
631             DO  k = nzb_2d(j,i) + 1, nzt
632                volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k)
633             ENDDO
634             !$OMP END CRITICAL
635          ENDIF
636          IF ( conserve_volume_flow  .AND.  bc_ns_cyc  .AND.  bc_lr_cyc  .AND. &
637               j == ny )  THEN
638             !$OMP CRITICAL
639             DO  k = nzb_2d(j,i) + 1, nzt
640                volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k)
641             ENDDO
642             !$OMP END CRITICAL
643          ENDIF
644
645       ENDDO
646    ENDDO
647    !$OMP END PARALLEL
648   
649!
650!-- Conserve the volume flow
651    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc  .AND.  bc_ns_cyc ) )  THEN
652
653#if defined( __parallel )   
654       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
655       CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 2, MPI_REAL, &
656                           MPI_SUM, comm2d, ierr ) 
657#else
658       volume_flow = volume_flow_l 
659#endif   
660
661       volume_flow_offset = ( volume_flow_initial - volume_flow ) / &
662                            volume_flow_area
663
664       !$OMP PARALLEL PRIVATE (i,j,k)
665       !$OMP DO
666       DO  i = nxl, nxr
667          DO  j = nys, nyn
668             DO  k = nzb_u_inner(j,i) + 1, nzt
669                u(k,j,i) = u(k,j,i) + volume_flow_offset(1)
670             ENDDO
671             DO k = nzb_v_inner(j,i) + 1, nzt
672                v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
673             ENDDO
674          ENDDO
675       ENDDO
676
677       !$OMP END PARALLEL
678
679    ENDIF
680
681!
682!-- Exchange of boundaries for the velocities
683    CALL exchange_horiz( u, nbgp )
684    CALL exchange_horiz( v, nbgp )
685    CALL exchange_horiz( w, nbgp )
686
687!
688!-- Compute the divergence of the corrected velocity field,
689!-- a possible PE-sum is computed in flow_statistics
690    CALL cpu_log( log_point_s(1), 'divergence', 'start' )
691    sums_divnew_l = 0.0
692
693!
694!-- d must be reset to zero because it can contain nonzero values below the
695!-- topography
696    IF ( topography /= 'flat' )  d = 0.0
697
698    localsum  = 0.0
699    threadsum = 0.0
700
701    !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
702    !$OMP DO SCHEDULE( STATIC )
703#if defined( __ibm )
704    DO  i = nxl, nxr
705       DO  j = nys, nyn
706          DO  k = nzb_s_inner(j,i)+1, nzt
707             d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
708                        ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
709                        ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
710          ENDDO
711          DO  k = nzb+1, nzt
712             threadsum = threadsum + ABS( d(k,j,i) )
713          ENDDO
714       ENDDO
715    ENDDO
716#else
717    DO  i = nxl, nxr
718       DO  j = nys, nyn
719          DO  k = nzb_s_inner(j,i)+1, nzt
720             d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
721                        ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
722                        ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
723             threadsum = threadsum + ABS( d(k,j,i) )
724          ENDDO
725       ENDDO
726    ENDDO
727#endif
728
729    localsum = localsum + threadsum
730    !$OMP END PARALLEL
731
732!
733!-- For completeness, set the divergence sum of all statistic regions to those
734!-- of the total domain
735    sums_divnew_l(0:statistic_regions) = localsum
736
737    CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
738
739    CALL cpu_log( log_point(8), 'pres', 'stop' )
740   
741
742
743 END SUBROUTINE pres
Note: See TracBrowser for help on using the repository browser.