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

Last change on this file since 4180 was 4180, checked in by scharf, 2 years ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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