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

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

code has been put under the GNU General Public License (v3)

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