source: palm/trunk/SOURCE/pres.f90

Last change on this file was 4855, checked in by raasch, 8 months ago

bugfix: mean w removal not applied to ghost points of the total domain in case of non-cyclic setups (pres), bugfix for correct identification of indices of extreme values in case of non-cyclic boundary conditions

  • Property svn:keywords set to Id
File size: 32.5 KB
Line 
1!> @file pres.f90
2!--------------------------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
6! Public License as published by the Free Software Foundation, either version 3 of the License, or
7! (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11! Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with PALM. If not, see
14! <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2021 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------------------------!
18!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: pres.f90 4855 2021-01-25 12:30:54Z banzhafs $
27! bugfix: mean w removal not applied to ghost points of the total domain in case of non-cyclic
28! setups
29!
30! 4828 2021-01-05 11:21:41Z Giersch
31! Bugfix in OpenMP directives - intel compiler do not allow reduction operations on array elements
32!
33! 4719 2020-10-01 11:28:33Z pavelkrc
34! Bugfix of previous commit (div_old discrepancy)
35!
36! 4717 2020-09-30 22:27:40Z pavelkrc
37! Fixes and optimizations of OpenMP parallelization, formatting of OpenMP
38! directives (J. Resler)
39!
40! 4651 2020-08-27 07:17:45Z raasch
41! preprocessor branch for ibm removed
42!
43! 4649 2020-08-25 12:11:17Z raasch
44! File re-formatted to follow the PALM coding standard
45!
46! 4457 2020-03-11 14:20:43Z raasch
47! use statement for exchange horiz added
48!
49! 4360 2020-01-07 11:25:50Z suehring
50! Introduction of wall_flags_total_0, which currently sets bits based on static
51! topography information used in wall_flags_static_0
52!
53! 4329 2019-12-10 15:46:36Z motisi
54! Renamed wall_flags_0 to wall_flags_static_0
55!
56! 4182 2019-08-22 15:20:23Z scharf
57! Corrected "Former revisions" section
58!
59! 4015 2019-06-05 13:25:35Z raasch
60! variable child_domain_nvn eliminated
61!
62! 3849 2019-04-01 16:35:16Z knoop
63! OpenACC port for SPEC
64!
65! Revision 1.1  1997/07/24 11:24:44  raasch
66! Initial revision
67!
68!
69!--------------------------------------------------------------------------------------------------!
70! Description:
71! ------------
72!> Compute the divergence of the provisional velocity field. Solve the Poisson equation for the
73!> perturbation pressure. Compute the final velocities using this perturbation pressure. Compute the
74!> remaining divergence.
75!--------------------------------------------------------------------------------------------------!
76 SUBROUTINE pres
77
78
79    USE arrays_3d,                                                                                 &
80        ONLY:  d,                                                                                  &
81               ddzu,                                                                               &
82               ddzu_pres,                                                                          &
83               ddzw,                                                                               &
84               dzw,                                                                                &
85               p,                                                                                  &
86               p_loc,                                                                              &
87               rho_air,                                                                            &
88               rho_air_zw,                                                                         &
89               tend,                                                                               &
90               u,                                                                                  &
91               v,                                                                                  &
92               w
93
94    USE control_parameters,                                                                        &
95        ONLY:  bc_lr_cyc,                                                                          &
96               bc_ns_cyc,                                                                          &
97               bc_radiation_l,                                                                     &
98               bc_radiation_n,                                                                     &
99               bc_radiation_r,                                                                     &
100               bc_radiation_s,                                                                     &
101               child_domain,                                                                       &
102               conserve_volume_flow,                                                               &
103               coupling_mode,                                                                      &
104               dt_3d,                                                                              &
105               gathered_size,                                                                      &
106               ibc_p_b,                                                                            &
107               ibc_p_t,                                                                            &
108               intermediate_timestep_count,                                                        &
109               intermediate_timestep_count_max,                                                    &
110               mg_switch_to_pe0_level,                                                             &
111               nesting_offline,                                                                    &
112               psolver,                                                                            &
113               subdomain_size,                                                                     &
114               topography,                                                                         &
115               volume_flow,                                                                        &
116               volume_flow_area,                                                                   &
117               volume_flow_initial
118
119    USE cpulog,                                                                                    &
120        ONLY:  cpu_log,                                                                            &
121               log_point,                                                                          &
122               log_point_s
123
124    USE exchange_horiz_mod,                                                                        &
125        ONLY:  exchange_horiz
126
127    USE grid_variables,                                                                            &
128        ONLY:  ddx,                                                                                &
129               ddy
130
131    USE indices,                                                                                   &
132        ONLY:  nbgp,                                                                               &
133               ngp_2dh_outer,                                                                      &
134               nx,                                                                                 &
135               nxl,                                                                                &
136               nxlg,                                                                               &
137               nxl_mg,                                                                             &
138               nxr,                                                                                &
139               nxrg,                                                                               &
140               nxr_mg,                                                                             &
141               ny,                                                                                 &
142               nys,                                                                                &
143               nysg,                                                                               &
144               nys_mg,                                                                             &
145               nyn,                                                                                &
146               nyng,                                                                               &
147               nyn_mg,                                                                             &
148               nzb,                                                                                &
149               nzt,                                                                                &
150               nzt_mg,                                                                             &
151               wall_flags_total_0
152
153    USE kinds
154
155    USE pegrid
156
157    USE pmc_interface,                                                                             &
158        ONLY:  nesting_mode
159
160    USE poisfft_mod,                                                                               &
161        ONLY:  poisfft
162
163    USE poismg_mod
164
165    USE poismg_noopt_mod
166
167    USE statistics,                                                                                &
168        ONLY:  statistic_regions,                                                                  &
169               sums_divnew_l,                                                                      &
170               sums_divold_l,                                                                      &
171               weight_pres,                                                                        &
172               weight_substep
173
174    USE surface_mod,                                                                               &
175        ONLY :  bc_h
176
177    IMPLICIT NONE
178
179    INTEGER(iwp) ::  i  !<
180    INTEGER(iwp) ::  j  !<
181    INTEGER(iwp) ::  k  !<
182    INTEGER(iwp) ::  m  !<
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
197    CALL cpu_log( log_point(8), 'pres', 'start' )
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 cannot be used in the
224!--    iterative solver. Therefore, its initial value is stored on p_loc, which is then iteratively
225!--    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 cannot be used in the
240!--    iterative solver. Therefore, its initial value is stored on p_loc, which is then iteratively
241!--    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 boundary conditions
248!-- WARNING: so far, this conservation does not work at the left/south boundary if the topography at
249!--          the inflow differs from that at the outflow! For this case, volume_flow_area needs
250!--          adjustment!
251!
252!-- Left/right
253    IF ( conserve_volume_flow  .AND.  ( bc_radiation_l .OR. bc_radiation_r ) )  THEN
254
255       volume_flow(1)   = 0.0_wp
256       volume_flow_l(1) = 0.0_wp
257
258       IF ( bc_radiation_l )  THEN
259          i = 0
260       ELSEIF ( bc_radiation_r )  THEN
261          i = nx+1
262       ENDIF
263
264       DO  j = nys, nyn
265!
266!--       Sum up the volume flow through the south/north boundary
267          DO  k = nzb+1, nzt
268             volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k)                               &
269                                * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
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, MPI_SUM, comm1dy, ierr )
276#else
277       volume_flow = volume_flow_l
278#endif
279       volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) ) / volume_flow_area(1)
280
281       DO  j = nysg, nyng
282          DO  k = nzb+1, nzt
283             u(k,j,i) = u(k,j,i) + volume_flow_offset(1)                                           &
284                        * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
285          ENDDO
286       ENDDO
287
288    ENDIF
289
290!
291!-- South/north
292    IF ( conserve_volume_flow  .AND.  ( bc_radiation_n .OR. bc_radiation_s ) )  THEN
293
294       volume_flow(2)   = 0.0_wp
295       volume_flow_l(2) = 0.0_wp
296
297       IF ( bc_radiation_s )  THEN
298          j = 0
299       ELSEIF ( bc_radiation_n )  THEN
300          j = ny+1
301       ENDIF
302
303       DO  i = nxl, nxr
304!
305!--       Sum up the volume flow through the south/north boundary
306          DO  k = nzb+1, nzt
307             volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k)                               &
308                                * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
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, MPI_SUM, comm1dx, ierr )
315#else
316       volume_flow = volume_flow_l
317#endif
318       volume_flow_offset(2) = ( volume_flow_initial(2) - volume_flow(2) ) / volume_flow_area(2)
319
320       DO  i = nxlg, nxrg
321          DO  k = nzb+1, nzt
322             v(k,j,i) = v(k,j,i) + volume_flow_offset(2)                                           &
323                        * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
324          ENDDO
325       ENDDO
326
327    ENDIF
328
329!
330!-- Remove mean vertical velocity in case that Neumann conditions are used both at bottom and top
331!-- boundary. With Neumann conditions at both vertical boundaries, the solver cannot remove mean
332!-- vertical velocities. They should be removed, because incompressibility requires that the
333!-- vertical gradient of vertical velocity is zero. Since w=0 at the solid surface, it must be zero
334!-- everywhere.
335!-- This must not be done in case of a 3d-nesting child domain, because a mean vertical velocity
336!-- can physically exist in such a domain.
337!-- Also in case of offline nesting, mean vertical velocities may exist (and must not be removed),
338!-- caused by horizontal divergence/convergence of the large scale flow that is prescribed at the
339!-- side boundaries.
340!-- The removal cannot be done before the first initial time step because ngp_2dh_outer is not yet
341!-- known then.
342    IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1  .AND.  .NOT. nesting_offline                           &
343         .AND. .NOT. ( child_domain .AND. nesting_mode /= 'vertical' )                             &
344         .AND. intermediate_timestep_count /= 0 )  THEN
345       w_l = 0.0_wp;  w_l_l = 0.0_wp
346       DO  i = nxl, nxr
347          DO  j = nys, nyn
348             DO  k = nzb+1, nzt
349                w_l_l(k) = w_l_l(k) + w(k,j,i)                                                     &
350                           * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
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, comm2d, ierr )
357#else
358       w_l = w_l_l
359#endif
360       DO  k = 1, nzt
361          w_l(k) = w_l(k) / ngp_2dh_outer(k,0)
362       ENDDO
363       DO  i = nxl, nxr
364          DO  j = nys, nyn
365             DO  k = nzb+1, nzt
366                w(k,j,i) = w(k,j,i) - w_l(k)                                                       &
367                           * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
368             ENDDO
369          ENDDO
370       ENDDO
371!
372!--    Instead of running the above loop over ghost points, they are set via exchange_horiz,
373!--    in order to correctly consider non-cyclic boundary conditions, where ghost boundaries
374!--    of the total domain must not be set. Otherwise w may continuously increase/decrease
375!--    at these points.
376       CALL exchange_horiz( w, nbgp )
377    ENDIF
378
379!
380!-- Compute the divergence of the provisional velocity field.
381    CALL cpu_log( log_point_s(1), 'divergence', 'start' )
382
383    IF ( psolver(1:9) == 'multigrid' )  THEN
384       !$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
385       DO  i = nxl-1, nxr+1
386          DO  j = nys-1, nyn+1
387             DO  k = nzb, nzt+1
388                d(k,j,i) = 0.0_wp
389             ENDDO
390          ENDDO
391       ENDDO
392    ELSE
393       !$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
394       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
395       !$ACC PRESENT(d)
396       DO  i = nxl, nxr
397          DO  j = nys, nyn
398             DO  k = nzb+1, nzt
399                d(k,j,i) = 0.0_wp
400             ENDDO
401          ENDDO
402       ENDDO
403    ENDIF
404
405    localsum  = 0.0_wp
406    threadsum = 0.0_wp
407
408    !$OMP PARALLEL PRIVATE (i,j,k)
409    !$OMP DO SCHEDULE( STATIC )
410    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
411    !$ACC PRESENT(u, v, w, rho_air, rho_air_zw, ddzw, wall_flags_total_0) &
412    !$ACC PRESENT(d)
413    DO  i = nxl, nxr
414       DO  j = nys, nyn
415          DO  k = 1, nzt
416             d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx +                           &
417                          ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy +                           &
418                          ( w(k,j,i)   * rho_air_zw(k) - w(k-1,j,i) * rho_air_zw(k-1) )            &
419                          * ddzw(k) ) * ddt_3d * d_weight_pres                                     &
420                        * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
421          ENDDO
422       ENDDO
423    ENDDO
424    !$OMP END PARALLEL
425
426!
427!-- Compute possible PE-sum of divergences for flow_statistics. Carry out computation only at last
428!-- Runge-Kutta substep.
429    IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.                      &
430         intermediate_timestep_count == 0 )  THEN
431       !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
432       !$OMP DO SCHEDULE( STATIC )
433       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
434       !$ACC REDUCTION(+:threadsum) COPY(threadsum) &
435       !$ACC PRESENT(d)
436       DO  i = nxl, nxr
437          DO  j = nys, nyn
438             DO  k = nzb+1, nzt
439                threadsum = threadsum + ABS( d(k,j,i) )
440             ENDDO
441          ENDDO
442       ENDDO
443       localsum = localsum + threadsum * dt_3d * weight_pres_l
444       !$OMP END PARALLEL
445    ENDIF
446
447!
448!-- For completeness, set the divergence sum of all statistic regions to those of the total domain
449    IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.                      &
450         intermediate_timestep_count == 0 )  THEN
451       sums_divold_l(0:statistic_regions) = localsum
452    ENDIF
453
454    CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
455
456!
457!-- Compute the pressure perturbation solving the Poisson equation
458    IF ( psolver == 'poisfft' )  THEN
459
460!
461!--    Solve Poisson equation via FFT and solution of tridiagonal matrices
462       CALL poisfft( d )
463
464!
465!--    Store computed perturbation pressure and set boundary condition in z-direction
466       !$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
467       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
468       !$ACC PRESENT(d, tend)
469       DO  i = nxl, nxr
470          DO  j = nys, nyn
471             DO  k = nzb+1, nzt
472                tend(k,j,i) = d(k,j,i)
473             ENDDO
474          ENDDO
475       ENDDO
476
477!
478!--    Bottom boundary:
479!--    This condition is only required for internal output. The pressure gradient
480!--    (dp(nzb+1)-dp(nzb))/dz is not used anywhere else.
481       IF ( ibc_p_b == 1 )  THEN
482!
483!--       Neumann (dp/dz = 0). Using surface data type, first for non-natural surfaces, then for
484!--       natural and urban surfaces
485!--       Upward facing
486          !$OMP PARALLEL DO PRIVATE( m, i, j, k ) SCHEDULE( STATIC )
487          !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
488          !$ACC PRESENT(bc_h, tend)
489          DO  m = 1, bc_h(0)%ns
490             i = bc_h(0)%i(m)
491             j = bc_h(0)%j(m)
492             k = bc_h(0)%k(m)
493             tend(k-1,j,i) = tend(k,j,i)
494          ENDDO
495!
496!--       Downward facing
497          !$OMP PARALLEL DO PRIVATE( m, i, j, k ) SCHEDULE( STATIC )
498          !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
499          !$ACC PRESENT(bc_h, tend)
500          DO  m = 1, bc_h(1)%ns
501             i = bc_h(1)%i(m)
502             j = bc_h(1)%j(m)
503             k = bc_h(1)%k(m)
504             tend(k+1,j,i) = tend(k,j,i)
505          ENDDO
506
507       ELSE
508!
509!--       Dirichlet. Using surface data type, first for non-natural surfaces, then for natural and
510!--       urban surfaces
511!--       Upward facing
512          !$OMP PARALLEL DO PRIVATE( m, i, j, k ) SCHEDULE( STATIC )
513          DO  m = 1, bc_h(0)%ns
514             i = bc_h(0)%i(m)
515             j = bc_h(0)%j(m)
516             k = bc_h(0)%k(m)
517             tend(k-1,j,i) = 0.0_wp
518          ENDDO
519!
520!--       Downward facing
521          !$OMP PARALLEL DO PRIVATE( m, i, j, k ) SCHEDULE( STATIC )
522          DO  m = 1, bc_h(1)%ns
523             i = bc_h(1)%i(m)
524             j = bc_h(1)%j(m)
525             k = bc_h(1)%k(m)
526             tend(k+1,j,i) = 0.0_wp
527          ENDDO
528
529       ENDIF
530
531!
532!--    Top boundary
533       IF ( ibc_p_t == 1 )  THEN
534!
535!--       Neumann
536          !$OMP PARALLEL DO PRIVATE (i,j) SCHEDULE( STATIC )
537          DO  i = nxlg, nxrg
538             DO  j = nysg, nyng
539                tend(nzt+1,j,i) = tend(nzt,j,i)
540             ENDDO
541          ENDDO
542
543       ELSE
544!
545!--       Dirichlet
546          !$OMP PARALLEL DO PRIVATE (i,j) SCHEDULE( STATIC )
547          !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
548          !$ACC PRESENT(tend)
549          DO  i = nxlg, nxrg
550             DO  j = nysg, nyng
551                tend(nzt+1,j,i) = 0.0_wp
552             ENDDO
553          ENDDO
554
555       ENDIF
556
557!
558!--    Exchange boundaries for p
559       CALL exchange_horiz( tend, nbgp )
560
561    ELSEIF ( psolver == 'sor' )  THEN
562
563!
564!--    Solve Poisson equation for perturbation pressure using SOR-Red/Black scheme
565       CALL sor( d, ddzu_pres, ddzw, p_loc )
566       tend = p_loc
567
568    ELSEIF ( psolver(1:9) == 'multigrid' )  THEN
569
570!
571!--    Solve Poisson equation for perturbation pressure using Multigrid scheme, array tend is used
572!--    to store the residuals.
573
574!--    If the number of grid points of the gathered grid, which is collected on PE0, is larger than
575!--    the number of grid points of an PE, than array tend will be enlarged.
576       IF ( gathered_size > subdomain_size )  THEN
577          DEALLOCATE( tend )
578          ALLOCATE( tend(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg(                              &
579                    mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1,                    &
580                    nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg(mg_switch_to_pe0_level)+1) )
581       ENDIF
582
583       IF ( psolver == 'multigrid' )  THEN
584          CALL poismg( tend )
585       ELSE
586          CALL poismg_noopt( tend )
587       ENDIF
588
589       IF ( gathered_size > subdomain_size )  THEN
590          DEALLOCATE( tend )
591          ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
592       ENDIF
593
594!
595!--    Restore perturbation pressure on tend because this array is used further below to correct the
596!--    velocity fields
597       DO  i = nxl-1, nxr+1
598          DO  j = nys-1, nyn+1
599             DO  k = nzb, nzt+1
600                tend(k,j,i) = p_loc(k,j,i)
601             ENDDO
602          ENDDO
603       ENDDO
604
605    ENDIF
606
607!
608!-- Store perturbation pressure on array p, used for pressure data output.
609!-- Ghost layers are added in the output routines (except sor-method: see below)
610    IF ( intermediate_timestep_count <= 1 )  THEN
611       !$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
612       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
613       !$ACC PRESENT(p, tend)
614       DO  i = nxl-1, nxr+1
615          DO  j = nys-1, nyn+1
616             DO  k = nzb, nzt+1
617                p(k,j,i) = tend(k,j,i) * weight_substep_l
618             ENDDO
619          ENDDO
620       ENDDO
621
622    ELSEIF ( intermediate_timestep_count > 1 )  THEN
623       !$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
624       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
625       !$ACC PRESENT(p, tend)
626       DO  i = nxl-1, nxr+1
627          DO  j = nys-1, nyn+1
628             DO  k = nzb, nzt+1
629                p(k,j,i) = p(k,j,i) + tend(k,j,i) * weight_substep_l
630             ENDDO
631          ENDDO
632       ENDDO
633
634    ENDIF
635
636!
637!-- SOR-method needs ghost layers for the next timestep
638    IF ( psolver == 'sor' )  CALL exchange_horiz( p, nbgp )
639
640!
641!-- Correction of the provisional velocities with the current perturbation pressure just computed
642    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc .OR. bc_ns_cyc ) )  THEN
643       volume_flow_l(1) = 0.0_wp
644       volume_flow_l(2) = 0.0_wp
645    ENDIF
646!
647!-- Add pressure gradients to the velocity components. Note, the loops are running over the entire
648!-- model domain, even though, in case of non-cyclic boundaries u- and v-component are not
649!-- prognostic at i=0 and j=0, respectiveley. However, in case of Dirichlet boundary conditions for
650!-- the velocities, zero-gradient conditions for the pressure are set, so that no modification is
651!-- imposed at the boundaries.
652    !$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
653    !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k) &
654    !$ACC PRESENT(u, v, w, tend, ddzu, wall_flags_total_0)
655    DO  i = nxl, nxr
656       DO  j = nys, nyn
657
658          DO  k = nzb+1, nzt
659             w(k,j,i) = w(k,j,i) - dt_3d * ( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1)             &
660                        * weight_pres_l                                                            &
661                        * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
662          ENDDO
663
664          DO  k = nzb+1, nzt
665             u(k,j,i) = u(k,j,i) - dt_3d * ( tend(k,j,i) - tend(k,j,i-1) ) * ddx * weight_pres_l   &
666                        * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
667          ENDDO
668
669          DO  k = nzb+1, nzt
670             v(k,j,i) = v(k,j,i) - dt_3d * ( tend(k,j,i) - tend(k,j-1,i) ) * ddy * weight_pres_l   &
671                        * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
672          ENDDO
673
674       ENDDO
675    ENDDO
676
677!
678!-- The vertical velocity is not set to zero at nzt + 1 for nested domains. Instead it is set to the
679!-- values of nzt (see routine vnest_boundary_conds or pmci_interp_tril_t) BEFORE calling the
680!-- pressure solver. To avoid jumps while plotting profiles, w at the top has to be set to the
681!-- values in height nzt after above modifications. Hint: w level nzt+1 does not impact results.
682    IF ( child_domain  .OR.  coupling_mode == 'vnested_fine' )  THEN
683       w(nzt+1,:,:) = w(nzt,:,:)
684    ENDIF
685
686!
687!-- Sum up the volume flow through the right and north boundary
688    IF ( conserve_volume_flow  .AND.  bc_lr_cyc  .AND.  bc_ns_cyc  .AND.  nxr == nx )  THEN
689       threadsum = 0.0_wp
690!
691!--    Summation of the volume flow is done on threadsum rather than on volumen_flow itself.
692!--    This is because intel compiler, when compiled with openmp, do not allow reduction
693!--    operation on array elements.
694       !$OMP PARALLEL DO PRIVATE (j,k) REDUCTION (+:threadsum)
695       DO  j = nys, nyn
696          DO  k = nzb+1, nzt
697             threadsum = threadsum + u(k,j,nxr) * dzw(k)                             &
698                                * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,nxr), 1 ) )
699          ENDDO
700       ENDDO
701       volume_flow_l(1) = threadsum
702
703    ENDIF
704
705    IF ( conserve_volume_flow  .AND.  bc_ns_cyc  .AND.  bc_lr_cyc  .AND. nyn == ny )  THEN
706       threadsum = 0.0_wp
707       !$OMP PARALLEL DO PRIVATE (j,k) REDUCTION (+:threadsum)
708       DO  i = nxl, nxr
709          DO  k = nzb+1, nzt
710             threadsum = threadsum + v(k,nyn,i) * dzw(k)                             &
711                                * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,nyn,i), 2 ) )
712           ENDDO
713       ENDDO
714       volume_flow_l(2) = threadsum
715
716    ENDIF
717
718!
719!-- Conserve the volume flow
720    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc  .AND.  bc_ns_cyc ) )  THEN
721
722#if defined( __parallel )
723       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
724       CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 2, MPI_REAL, MPI_SUM, comm2d, ierr )
725#else
726       volume_flow = volume_flow_l
727#endif
728
729       volume_flow_offset(1:2) = ( volume_flow_initial(1:2) - volume_flow(1:2) )                   &
730                                 / volume_flow_area(1:2)
731
732       !$OMP PARALLEL DO PRIVATE (i,j,k) SCHEDULE( STATIC )
733       DO  i = nxl, nxr
734          DO  j = nys, nyn
735             DO  k = nzb+1, nzt
736                u(k,j,i) = u(k,j,i) + volume_flow_offset(1)                                        &
737                           * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
738             ENDDO
739             DO  k = nzb+1, nzt
740                v(k,j,i) = v(k,j,i) + volume_flow_offset(2)                                        &
741                           * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
742             ENDDO
743          ENDDO
744       ENDDO
745
746    ENDIF
747
748!
749!-- Exchange of boundaries for the velocities
750    CALL exchange_horiz( u, nbgp )
751    CALL exchange_horiz( v, nbgp )
752    CALL exchange_horiz( w, nbgp )
753
754!
755!-- Compute the divergence of the corrected velocity field.
756!-- A possible PE-sum is computed in flow_statistics. Carry out computation only at last
757!-- Runge-Kutta step.
758    IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.                      &
759         intermediate_timestep_count == 0 )  THEN
760       CALL cpu_log( log_point_s(1), 'divergence', 'start' )
761       sums_divnew_l = 0.0_wp
762
763!
764!--    d must be reset to zero because it can contain nonzero values below the topography
765       IF ( topography /= 'flat' )  d = 0.0_wp
766
767       localsum  = 0.0_wp
768       threadsum = 0.0_wp
769
770       !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
771       !$OMP DO SCHEDULE( STATIC )
772       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
773       !$ACC PRESENT(u, v, w, rho_air, rho_air_zw, ddzw, wall_flags_total_0) &
774       !$ACC PRESENT(d)
775       DO  i = nxl, nxr
776          DO  j = nys, nyn
777             DO  k = nzb+1, nzt
778                d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx +    &
779                             ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy +    &
780                             ( w(k,j,i)   * rho_air_zw(k) - w(k-1,j,i) * rho_air_zw(k-1) )         &
781                             * ddzw(k) )                                                           &
782                           * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
783             ENDDO
784          ENDDO
785       ENDDO
786!
787!--    Compute possible PE-sum of divergences for flow_statistics
788       !$OMP DO SCHEDULE( STATIC )
789       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
790       !$ACC REDUCTION(+:threadsum) COPY(threadsum) &
791       !$ACC PRESENT(d)
792       DO  i = nxl, nxr
793          DO  j = nys, nyn
794             DO  k = nzb+1, nzt
795                threadsum = threadsum + ABS( d(k,j,i) )
796             ENDDO
797          ENDDO
798       ENDDO
799
800       localsum = localsum + threadsum
801       !$OMP END PARALLEL
802
803!
804!--    For completeness, set the divergence sum of all statistic regions to those of the total
805!--    domain
806       sums_divnew_l(0:statistic_regions) = localsum
807
808       CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
809
810    ENDIF
811
812    CALL cpu_log( log_point(8), 'pres', 'stop' )
813
814 END SUBROUTINE pres
Note: See TracBrowser for help on using the repository browser.