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

Last change on this file since 2118 was 2118, checked in by raasch, 5 years ago

all OpenACC directives and related parts removed from the code

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