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

Last change on this file since 4752 was 4740, checked in by suehring, 4 years ago

Bugfix in OpenMP directives - intel compiler do not allow reduction operations on array elements

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