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

Last change on this file since 2186 was 2119, checked in by raasch, 8 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 25.7 KB
Line 
1!> @file pres.f90
2!------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! 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-2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: pres.f90 2119 2017-01-17 16:51:50Z raasch $
27!
28! 2118 2017-01-17 16:38:49Z raasch
29! OpenACC directives and related code removed
30!
31! 2073 2016-11-30 14:34:05Z raasch
32! openmp bugfix for calculation of new divergence
33!
34! 2037 2016-10-26 11:15:40Z knoop
35! Anelastic approximation implemented
36!
37! 2000 2016-08-20 18:09:15Z knoop
38! Forced header and separation lines into 80 columns
39!
40! 1932 2016-06-10 12:09:21Z suehring
41! Initial version of purely vertical nesting introduced.
42!
43! 1931 2016-06-10 12:06:59Z suehring
44! Rename multigrid into multigrid_noopt and multigrid_fast into multigrid
45!
46! 1929 2016-06-09 16:25:25Z suehring
47! Bugfix: weight_substep for initial call, replace by local variable
48!
49! 1918 2016-05-27 14:35:57Z raasch
50! sum of divergence is also calculated when pres is called before the initial
51! first time step,
52! stearing is modified by using intermediate_timestep_count = 0 in order to
53! determine if pres is called before the first initial timestep,
54! bugfix: in case of Neumann conditions for pressure both at bottom and top,
55!         mean vertical velocity is also removed for the first time step
56! bugfix for calculating divergences
57!
58! 1908 2016-05-25 17:22:32Z suehring
59! New divergence for output into RUN_CONTROL file is calculated only at last
60! Runge-Kutta step
61!
62! 1845 2016-04-08 08:29:13Z raasch
63! nzb_2d replace by nzb_u|v_inner
64!
65! 1799 2016-04-05 08:35:55Z gronemeier
66! Bugfix: excluded third dimension from horizontal volume flow calculation
67!
68! 1762 2016-02-25 12:31:13Z hellstea
69! Introduction of nested domain feature
70!
71! 1682 2015-10-07 23:56:08Z knoop
72! Code annotations made doxygen readable
73!
74! 1575 2015-03-27 09:56:27Z raasch
75! poismg_fast + respective module added, adjustments for psolver-queries
76!
77! 1342 2014-03-26 17:04:47Z kanani
78! REAL constants defined as wp-kind
79!
80! 1320 2014-03-20 08:40:49Z raasch
81! ONLY-attribute added to USE-statements,
82! kind-parameters added to all INTEGER and REAL declaration statements,
83! kinds are defined in new module kinds,
84! old module precision_kind is removed,
85! revision history before 2012 removed,
86! comment fields (!:) to be used for variable explanations added to
87! all variable declaration statements
88!
89! 1318 2014-03-17 13:35:16Z raasch
90! module interfaces removed
91!
92! 1306 2014-03-13 14:30:59Z raasch
93! second argument removed from routine poisfft
94!
95! 1257 2013-11-08 15:18:40Z raasch
96! openacc loop and loop vector clauses removed, independent clauses added,
97! end parallel replaced by end parallel loop
98!
99! 1221 2013-09-10 08:59:13Z raasch
100! openACC porting of reduction operations, loops for calculating d are
101! using the rflags_s_inner multiply flag instead of the nzb_s_inner loop index
102!
103! 1212 2013-08-15 08:46:27Z raasch
104! call of poisfft_hybrid removed
105!
106! 1117 2013-03-27 11:15:36Z suehring
107! Bugfix in OpenMP parallelization.
108!
109! 1113 2013-03-10 02:48:14Z raasch
110! GPU-porting of several loops, some loops rearranged
111!
112! 1111 2013-03-08 23:54:10Z
113! openACC statements added,
114! ibc_p_b = 2 removed
115!
116! 1092 2013-02-02 11:24:22Z raasch
117! unused variables removed
118!
119! 1036 2012-10-22 13:43:42Z raasch
120! code put under GPL (PALM 3.9)
121!
122! 1003 2012-09-14 14:35:53Z raasch
123! adjustment of array tend for cases with unequal subdomain sizes removed
124!
125! Revision 1.1  1997/07/24 11:24:44  raasch
126! Initial revision
127!
128!
129! Description:
130! ------------
131!> Compute the divergence of the provisional velocity field. Solve the Poisson
132!> equation for the perturbation pressure. Compute the final velocities using
133!> this perturbation pressure. Compute the remaining divergence.
134!------------------------------------------------------------------------------!
135 SUBROUTINE pres
136 
137
138    USE arrays_3d,                                                             &
139        ONLY:  d, ddzu, ddzu_pres, ddzw, dzw, p, p_loc, rho_air, rho_air_zw,   &
140               tend, u, v, w
141
142    USE control_parameters,                                                    &
143        ONLY:  bc_lr_cyc, bc_ns_cyc, conserve_volume_flow, dt_3d,              &
144               gathered_size, ibc_p_b, ibc_p_t, intermediate_timestep_count,   &
145               intermediate_timestep_count_max, mg_switch_to_pe0_level,        &
146               nest_domain, outflow_l, outflow_n, outflow_r,                   &
147               outflow_s, psolver, subdomain_size, topography, volume_flow,    &
148               volume_flow_area, volume_flow_initial
149
150    USE cpulog,                                                                &
151        ONLY:  cpu_log, log_point, log_point_s
152
153    USE grid_variables,                                                        &
154        ONLY:  ddx, ddy
155
156    USE indices,                                                               &
157        ONLY:  nbgp, ngp_2dh_outer, nx, nxl, nxlg, nxl_mg, nxr, nxrg, nxr_mg,  &
158               ny, nys, nysg, nys_mg, nyn, nyng, nyn_mg, nzb, nzb_s_inner,     &
159               nzb_u_inner, nzb_v_inner, nzb_w_inner, nzt, nzt_mg,             &
160               rflags_s_inner
161
162    USE kinds
163
164    USE pegrid
165   
166    USE pmc_interface,                                                         &
167        ONLY:  nesting_mode 
168
169    USE poisfft_mod,                                                           &
170        ONLY:  poisfft
171
172    USE poismg_mod
173
174    USE statistics,                                                            &
175        ONLY:  statistic_regions, sums_divnew_l, sums_divold_l, weight_pres,   &
176               weight_substep
177
178    IMPLICIT NONE
179
180    INTEGER(iwp) ::  i              !<
181    INTEGER(iwp) ::  j              !<
182    INTEGER(iwp) ::  k              !<
183
184    REAL(wp)     ::  ddt_3d         !<
185    REAL(wp)     ::  d_weight_pres  !<
186    REAL(wp)     ::  localsum       !<
187    REAL(wp)     ::  threadsum      !<
188    REAL(wp)     ::  weight_pres_l  !<
189    REAL(wp)     ::  weight_substep_l !<
190
191    REAL(wp), DIMENSION(1:3)   ::  volume_flow_l       !<
192    REAL(wp), DIMENSION(1:3)   ::  volume_flow_offset  !<
193    REAL(wp), DIMENSION(1:nzt) ::  w_l                 !<
194    REAL(wp), DIMENSION(1:nzt) ::  w_l_l               !<
195
196    LOGICAL :: nest_domain_nvn      !<
197
198
199    CALL cpu_log( log_point(8), 'pres', 'start' )
200
201
202!
203!-- Calculate quantities to be used locally
204    ddt_3d = 1.0_wp / dt_3d
205    IF ( intermediate_timestep_count == 0 )  THEN
206!
207!--    If pres is called before initial time step
208       weight_pres_l    = 1.0_wp
209       d_weight_pres    = 1.0_wp
210       weight_substep_l = 1.0_wp
211    ELSE
212       weight_pres_l    = weight_pres(intermediate_timestep_count)
213       d_weight_pres    = 1.0_wp / weight_pres(intermediate_timestep_count)
214       weight_substep_l = weight_substep(intermediate_timestep_count)
215    ENDIF
216
217!
218!-- Multigrid method expects array d to have one ghost layer.
219!--
220    IF ( psolver(1:9) == 'multigrid' )  THEN
221     
222       DEALLOCATE( d )
223       ALLOCATE( d(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ) 
224
225!
226!--    Since p is later used to hold the weighted average of the substeps, it
227!--    cannot be used in the iterative solver. Therefore, its initial value is
228!--    stored on p_loc, which is then iteratively advanced in every substep.
229       IF ( intermediate_timestep_count <= 1 )  THEN
230          DO  i = nxl-1, nxr+1
231             DO  j = nys-1, nyn+1
232                DO  k = nzb, nzt+1
233                   p_loc(k,j,i) = p(k,j,i)
234                ENDDO
235             ENDDO
236          ENDDO
237       ENDIF
238       
239    ELSEIF ( psolver == 'sor'  .AND.  intermediate_timestep_count <= 1 )  THEN
240
241!
242!--    Since p is later used to hold the weighted average of the substeps, it
243!--    cannot be used in the iterative solver. Therefore, its initial value is
244!--    stored on p_loc, which is then iteratively advanced in every substep.
245       p_loc = p
246
247    ENDIF
248
249!
250!-- Conserve the volume flow at the outflow in case of non-cyclic lateral
251!-- boundary conditions
252!-- WARNING: so far, this conservation does not work at the left/south
253!--          boundary if the topography at the inflow differs from that at the
254!--          outflow! For this case, volume_flow_area needs adjustment!
255!
256!-- Left/right
257    IF ( conserve_volume_flow  .AND.  ( outflow_l .OR. outflow_r ) )  THEN
258
259       volume_flow(1)   = 0.0_wp
260       volume_flow_l(1) = 0.0_wp
261
262       IF ( outflow_l )  THEN
263          i = 0
264       ELSEIF ( outflow_r )  THEN
265          i = nx+1
266       ENDIF
267
268       DO  j = nys, nyn
269!
270!--       Sum up the volume flow through the south/north boundary
271          DO  k = nzb_u_inner(j,i)+1, nzt
272             volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k)
273          ENDDO
274       ENDDO
275
276#if defined( __parallel )   
277       IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
278       CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, &
279                           MPI_SUM, comm1dy, ierr )   
280#else
281       volume_flow = volume_flow_l 
282#endif
283       volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) ) &
284                               / volume_flow_area(1)
285
286       DO  j = nysg, nyng
287          DO  k = nzb_u_inner(j,i)+1, nzt
288             u(k,j,i) = u(k,j,i) + volume_flow_offset(1)
289          ENDDO
290       ENDDO
291
292    ENDIF
293
294!
295!-- South/north
296    IF ( conserve_volume_flow  .AND.  ( outflow_n .OR. outflow_s ) )  THEN
297
298       volume_flow(2)   = 0.0_wp
299       volume_flow_l(2) = 0.0_wp
300
301       IF ( outflow_s )  THEN
302          j = 0
303       ELSEIF ( outflow_n )  THEN
304          j = ny+1
305       ENDIF
306
307       DO  i = nxl, nxr
308!
309!--       Sum up the volume flow through the south/north boundary
310          DO  k = nzb_v_inner(j,i)+1, nzt
311             volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k)
312          ENDDO
313       ENDDO
314
315#if defined( __parallel )   
316       IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
317       CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL, &
318                           MPI_SUM, comm1dx, ierr )   
319#else
320       volume_flow = volume_flow_l 
321#endif
322       volume_flow_offset(2) = ( volume_flow_initial(2) - volume_flow(2) )    &
323                               / volume_flow_area(2)
324
325       DO  i = nxlg, nxrg
326          DO  k = nzb_v_inner(j,i)+1, nzt
327             v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
328          ENDDO
329       ENDDO
330
331    ENDIF
332
333!
334!-- Remove mean vertical velocity in case that Neumann conditions are
335!-- used both at bottom and top boundary, and if not a nested domain in a
336!-- normal nesting run. In case of vertical nesting, this must be done.
337!-- Therefore an auxiliary logical variable nest_domain_nvn is used here, and
338!-- nvn stands for non-vertical nesting.
339!-- This cannot be done before the first initial time step because ngp_2dh_outer
340!-- is not yet known then.
341    nest_domain_nvn = nest_domain
342    IF ( nest_domain .AND. nesting_mode == 'vertical' )  THEN
343       nest_domain_nvn = .FALSE.
344    ENDIF
345
346    IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1  .AND.                               &
347         .NOT. nest_domain_nvn  .AND. intermediate_timestep_count /= 0 )        &
348    THEN
349       w_l = 0.0_wp;  w_l_l = 0.0_wp
350       DO  i = nxl, nxr
351          DO  j = nys, nyn
352             DO  k = nzb_w_inner(j,i)+1, nzt
353                w_l_l(k) = w_l_l(k) + w(k,j,i)
354             ENDDO
355          ENDDO
356       ENDDO
357#if defined( __parallel )   
358       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
359       CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, &
360                           comm2d, ierr )
361#else
362       w_l = w_l_l
363#endif
364       DO  k = 1, nzt
365          w_l(k) = w_l(k) / ngp_2dh_outer(k,0)
366       ENDDO
367       DO  i = nxlg, nxrg
368          DO  j = nysg, nyng
369             DO  k = nzb_w_inner(j,i)+1, nzt
370                w(k,j,i) = w(k,j,i) - w_l(k)
371             ENDDO
372          ENDDO
373       ENDDO
374    ENDIF
375
376!
377!-- Compute the divergence of the provisional velocity field.
378    CALL cpu_log( log_point_s(1), 'divergence', 'start' )
379
380    IF ( psolver(1:9) == 'multigrid' )  THEN
381       !$OMP PARALLEL DO SCHEDULE( STATIC )
382       DO  i = nxl-1, nxr+1
383          DO  j = nys-1, nyn+1
384             DO  k = nzb, nzt+1
385                d(k,j,i) = 0.0_wp
386             ENDDO
387          ENDDO
388       ENDDO
389    ELSE
390       !$OMP PARALLEL DO SCHEDULE( STATIC )
391       DO  i = nxl, nxr
392          DO  j = nys, nyn
393             DO  k = nzb+1, nzt
394                d(k,j,i) = 0.0_wp
395             ENDDO
396          ENDDO
397       ENDDO
398    ENDIF
399
400    localsum  = 0.0_wp
401    threadsum = 0.0_wp
402
403#if defined( __ibm )
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_s_inner(j,i)+1, nzt
409             d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx +       &
410                          ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy +       &
411                          ( w(k,j,i)   * rho_air_zw(k) -                       &
412                            w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k)           &
413                        ) * ddt_3d * d_weight_pres 
414          ENDDO
415!
416!--       Compute possible PE-sum of divergences for flow_statistics
417          DO  k = nzb_s_inner(j,i)+1, nzt
418             threadsum = threadsum + ABS( d(k,j,i) )
419          ENDDO
420
421       ENDDO
422    ENDDO
423
424    IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.  &
425         intermediate_timestep_count == 0 )  THEN
426       localsum = localsum + threadsum * dt_3d * weight_pres_l
427    ENDIF
428    !$OMP END PARALLEL
429#else
430
431    !$OMP PARALLEL PRIVATE (i,j,k)
432    !$OMP DO SCHEDULE( STATIC )
433    DO  i = nxl, nxr
434       DO  j = nys, nyn
435          DO  k = 1, nzt
436             d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx +       &
437                          ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy +       &
438                          ( w(k,j,i)   * rho_air_zw(k) -                       &
439                            w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k)           &
440                        ) * ddt_3d * d_weight_pres * rflags_s_inner(k,j,i)
441          ENDDO
442       ENDDO
443    ENDDO
444    !$OMP END PARALLEL
445
446!
447!-- Compute possible PE-sum of divergences for flow_statistics. Carry out
448!-- computation only at last Runge-Kutta substep.
449    IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.  &
450         intermediate_timestep_count == 0 )  THEN
451       !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
452       !$OMP DO SCHEDULE( STATIC )
453       DO  i = nxl, nxr
454          DO  j = nys, nyn
455             DO  k = nzb+1, nzt
456                threadsum = threadsum + ABS( d(k,j,i) )
457             ENDDO
458          ENDDO
459       ENDDO
460       localsum = localsum + threadsum * dt_3d * weight_pres_l
461       !$OMP END PARALLEL
462    ENDIF
463#endif
464
465!
466!-- For completeness, set the divergence sum of all statistic regions to those
467!-- of the total domain
468    IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.  &
469         intermediate_timestep_count == 0 )  THEN
470       sums_divold_l(0:statistic_regions) = localsum
471    ENDIF
472
473    CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
474
475!
476!-- Compute the pressure perturbation solving the Poisson equation
477    IF ( psolver == 'poisfft' )  THEN
478
479!
480!--    Solve Poisson equation via FFT and solution of tridiagonal matrices
481       CALL poisfft( d )
482
483!
484!--    Store computed perturbation pressure and set boundary condition in
485!--    z-direction
486       !$OMP PARALLEL DO
487       DO  i = nxl, nxr
488          DO  j = nys, nyn
489             DO  k = nzb+1, nzt
490                tend(k,j,i) = d(k,j,i)
491             ENDDO
492          ENDDO
493       ENDDO
494
495!
496!--    Bottom boundary:
497!--    This condition is only required for internal output. The pressure
498!--    gradient (dp(nzb+1)-dp(nzb))/dz is not used anywhere else.
499       IF ( ibc_p_b == 1 )  THEN
500!
501!--       Neumann (dp/dz = 0)
502          !$OMP PARALLEL DO
503          DO  i = nxlg, nxrg
504             DO  j = nysg, nyng
505                tend(nzb_s_inner(j,i),j,i) = tend(nzb_s_inner(j,i)+1,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(nzb_s_inner(j,i),j,i) = 0.0_wp
516             ENDDO
517          ENDDO
518
519       ENDIF
520
521!
522!--    Top boundary
523       IF ( ibc_p_t == 1 )  THEN
524!
525!--       Neumann
526          !$OMP PARALLEL DO
527          DO  i = nxlg, nxrg
528             DO  j = nysg, nyng
529                tend(nzt+1,j,i) = tend(nzt,j,i)
530             ENDDO
531          ENDDO
532
533       ELSE
534!
535!--       Dirichlet
536          !$OMP PARALLEL DO
537          DO  i = nxlg, nxrg
538             DO  j = nysg, nyng
539                tend(nzt+1,j,i) = 0.0_wp
540             ENDDO
541          ENDDO
542
543       ENDIF
544
545!
546!--    Exchange boundaries for p
547       CALL exchange_horiz( tend, nbgp )
548     
549    ELSEIF ( psolver == 'sor' )  THEN
550
551!
552!--    Solve Poisson equation for perturbation pressure using SOR-Red/Black
553!--    scheme
554       CALL sor( d, ddzu_pres, ddzw, p_loc )
555       tend = p_loc
556
557    ELSEIF ( psolver(1:9) == 'multigrid' )  THEN
558
559!
560!--    Solve Poisson equation for perturbation pressure using Multigrid scheme,
561!--    array tend is used to store the residuals, logical exchange_mg is used
562!--    to discern data exchange in multigrid ( 1 ghostpoint ) and normal grid
563!--    ( nbgp ghost points ).
564
565!--    If the number of grid points of the gathered grid, which is collected
566!--    on PE0, is larger than the number of grid points of an PE, than array
567!--    tend will be enlarged.
568       IF ( gathered_size > subdomain_size )  THEN
569          DEALLOCATE( tend )
570          ALLOCATE( tend(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg(          &
571                    mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1,&
572                    nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg(                   &
573                    mg_switch_to_pe0_level)+1) )
574       ENDIF
575
576       IF ( psolver == 'multigrid' )  THEN
577          CALL poismg( tend )
578       ELSE
579          CALL poismg_noopt( tend )
580       ENDIF
581
582       IF ( gathered_size > subdomain_size )  THEN
583          DEALLOCATE( tend )
584          ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
585       ENDIF
586
587!
588!--    Restore perturbation pressure on tend because this array is used
589!--    further below to correct the velocity fields
590       DO  i = nxl-1, nxr+1
591          DO  j = nys-1, nyn+1
592             DO  k = nzb, nzt+1
593                tend(k,j,i) = p_loc(k,j,i)
594             ENDDO
595          ENDDO
596       ENDDO
597
598    ENDIF
599
600!
601!-- Store perturbation pressure on array p, used for pressure data output.
602!-- Ghost layers are added in the output routines (except sor-method: see below)
603    IF ( intermediate_timestep_count <= 1 )  THEN
604       !$OMP PARALLEL PRIVATE (i,j,k)
605       !$OMP DO
606       DO  i = nxl-1, nxr+1
607          DO  j = nys-1, nyn+1
608             DO  k = nzb, nzt+1
609                p(k,j,i) = tend(k,j,i) * &
610                           weight_substep_l
611             ENDDO
612          ENDDO
613       ENDDO
614       !$OMP END PARALLEL
615
616    ELSEIF ( intermediate_timestep_count > 1 )  THEN
617       !$OMP PARALLEL PRIVATE (i,j,k)
618       !$OMP DO
619       DO  i = nxl-1, nxr+1
620          DO  j = nys-1, nyn+1
621             DO  k = nzb, nzt+1
622                p(k,j,i) = p(k,j,i) + tend(k,j,i) * &
623                           weight_substep_l
624             ENDDO
625          ENDDO
626       ENDDO
627       !$OMP END PARALLEL
628
629    ENDIF
630       
631!
632!-- SOR-method needs ghost layers for the next timestep
633    IF ( psolver == 'sor' )  CALL exchange_horiz( p, nbgp )
634
635!
636!-- Correction of the provisional velocities with the current perturbation
637!-- pressure just computed
638    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc .OR. bc_ns_cyc ) )  THEN
639       volume_flow_l(1) = 0.0_wp
640       volume_flow_l(2) = 0.0_wp
641    ENDIF
642
643    !$OMP PARALLEL PRIVATE (i,j,k)
644    !$OMP DO
645    DO  i = nxl, nxr   
646       DO  j = nys, nyn
647
648          DO  k = 1, nzt
649             IF ( k > nzb_w_inner(j,i) )  THEN
650                w(k,j,i) = w(k,j,i) - dt_3d *                                 &
651                           ( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1) *      &
652                           weight_pres_l
653             ENDIF
654          ENDDO
655
656          DO  k = 1, nzt
657             IF ( k > nzb_u_inner(j,i) )  THEN
658                u(k,j,i) = u(k,j,i) - dt_3d *                                 &
659                           ( tend(k,j,i) - tend(k,j,i-1) ) * ddx *            &
660                           weight_pres_l
661             ENDIF
662          ENDDO
663
664          DO  k = 1, nzt
665             IF ( k > nzb_v_inner(j,i) )  THEN
666                v(k,j,i) = v(k,j,i) - dt_3d *                                 &
667                           ( tend(k,j,i) - tend(k,j-1,i) ) * ddy *            &
668                           weight_pres_l
669             ENDIF
670          ENDDO                                                         
671
672       ENDDO
673    ENDDO
674    !$OMP END PARALLEL
675
676!
677!-- Sum up the volume flow through the right and north boundary
678    IF ( conserve_volume_flow  .AND.  bc_lr_cyc  .AND.  bc_ns_cyc  .AND.  &
679         nxr == nx )  THEN
680
681       !$OMP PARALLEL PRIVATE (j,k)
682       !$OMP DO
683       DO  j = nys, nyn
684          !$OMP CRITICAL
685          DO  k = nzb_u_inner(j,nx) + 1, nzt
686             volume_flow_l(1) = volume_flow_l(1) + u(k,j,nx) * dzw(k)
687          ENDDO
688          !$OMP END CRITICAL
689       ENDDO
690       !$OMP END PARALLEL
691
692    ENDIF
693
694    IF ( conserve_volume_flow  .AND.  bc_ns_cyc  .AND.  bc_lr_cyc  .AND.  &
695         nyn == ny )  THEN
696
697       !$OMP PARALLEL PRIVATE (i,k)
698       !$OMP DO
699       DO  i = nxl, nxr
700          !$OMP CRITICAL
701          DO  k = nzb_v_inner(ny,i) + 1, nzt
702             volume_flow_l(2) = volume_flow_l(2) + v(k,ny,i) * dzw(k)
703           ENDDO
704          !$OMP END CRITICAL
705       ENDDO
706       !$OMP END PARALLEL
707
708    ENDIF
709   
710!
711!-- Conserve the volume flow
712    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc  .AND.  bc_ns_cyc ) )  THEN
713
714#if defined( __parallel )   
715       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
716       CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 2, MPI_REAL, &
717                           MPI_SUM, comm2d, ierr ) 
718#else
719       volume_flow = volume_flow_l 
720#endif   
721
722       volume_flow_offset(1:2) = ( volume_flow_initial(1:2) - volume_flow(1:2) ) / &
723                            volume_flow_area(1:2)
724
725       !$OMP PARALLEL PRIVATE (i,j,k)
726       !$OMP DO
727       DO  i = nxl, nxr
728          DO  j = nys, nyn
729             DO  k = nzb_u_inner(j,i) + 1, nzt
730                u(k,j,i) = u(k,j,i) + volume_flow_offset(1)
731             ENDDO
732             DO  k = nzb_v_inner(j,i) + 1, nzt
733                v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
734             ENDDO
735          ENDDO
736       ENDDO
737
738       !$OMP END PARALLEL
739
740    ENDIF
741
742!
743!-- Exchange of boundaries for the velocities
744    CALL exchange_horiz( u, nbgp )
745    CALL exchange_horiz( v, nbgp )
746    CALL exchange_horiz( w, nbgp )
747
748!
749!-- Compute the divergence of the corrected velocity field,
750!-- A possible PE-sum is computed in flow_statistics. Carry out computation
751!-- only at last Runge-Kutta step.
752    IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.  &
753         intermediate_timestep_count == 0 )  THEN
754       CALL cpu_log( log_point_s(1), 'divergence', 'start' )
755       sums_divnew_l = 0.0_wp
756
757!
758!--    d must be reset to zero because it can contain nonzero values below the
759!--    topography
760       IF ( topography /= 'flat' )  d = 0.0_wp
761
762       localsum  = 0.0_wp
763       threadsum = 0.0_wp
764
765       !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
766#if defined( __ibm )
767       !$OMP DO SCHEDULE( STATIC )
768       DO  i = nxl, nxr
769          DO  j = nys, nyn
770             DO  k = nzb_s_inner(j,i)+1, nzt
771             d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx +         &
772                        ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy +         &
773                        ( w(k,j,i)   * rho_air_zw(k) -                         &
774                          w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k)
775             ENDDO
776             DO  k = nzb+1, nzt
777                threadsum = threadsum + ABS( d(k,j,i) )
778             ENDDO
779          ENDDO
780       ENDDO
781#else
782       !$OMP DO SCHEDULE( STATIC )
783       DO  i = nxl, nxr
784          DO  j = nys, nyn
785             DO  k = 1, nzt
786                d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx +    &
787                             ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy +    &
788                             ( w(k,j,i)   * rho_air_zw(k) -                    &
789                               w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k)        &
790                           ) * rflags_s_inner(k,j,i)
791             ENDDO
792          ENDDO
793       ENDDO
794!
795!--    Compute possible PE-sum of divergences for flow_statistics
796       !$OMP DO SCHEDULE( STATIC )
797       DO  i = nxl, nxr
798          DO  j = nys, nyn
799             DO  k = nzb+1, nzt
800                threadsum = threadsum + ABS( d(k,j,i) )
801             ENDDO
802          ENDDO
803       ENDDO
804#endif
805
806       localsum = localsum + threadsum
807       !$OMP END PARALLEL
808
809!
810!--    For completeness, set the divergence sum of all statistic regions to those
811!--    of the total domain
812       sums_divnew_l(0:statistic_regions) = localsum
813
814       CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
815
816    ENDIF
817
818    CALL cpu_log( log_point(8), 'pres', 'stop' )
819
820
821 END SUBROUTINE pres
Note: See TracBrowser for help on using the repository browser.