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

Last change on this file since 1094 was 1093, checked in by raasch, 11 years ago

last commit documented

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