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

Last change on this file since 4651 was 4651, checked in by raasch, 4 years ago

preprocessor branch for ibm removed

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