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

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

last commit documented

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