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

Last change on this file since 4729 was 4719, checked in by pavelkrc, 4 years ago

Bugfix of previous commit (div_old discrepancy)

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