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

Last change on this file since 1320 was 1320, checked in by raasch, 7 years ago

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

  • Property svn:keywords set to Id
File size: 23.8 KB
Line 
1 SUBROUTINE pres
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later 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-2014 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22! ONLY-attribute added to USE-statements,
23! kind-parameters added to all INTEGER and REAL declaration statements,
24! kinds are defined in new module kinds,
25! old module precision_kind is removed,
26! revision history before 2012 removed,
27! comment fields (!:) to be used for variable explanations added to
28! all variable declaration statements
29!
30! Former revisions:
31! -----------------
32! $Id: pres.f90 1320 2014-03-20 08:40:49Z raasch $
33!
34! 1318 2014-03-17 13:35:16Z raasch
35! module interfaces removed
36!
37! 1306 2014-03-13 14:30:59Z raasch
38! second argument removed from routine poisfft
39!
40! 1257 2013-11-08 15:18:40Z raasch
41! openacc loop and loop vector clauses removed, independent clauses added,
42! end parallel replaced by end parallel loop
43!
44! 1221 2013-09-10 08:59:13Z raasch
45! openACC porting of reduction operations, loops for calculating d are
46! using the rflags_s_inner multiply flag instead of the nzb_s_inner loop index
47!
48! 1212 2013-08-15 08:46:27Z raasch
49! call of poisfft_hybrid removed
50!
51! 1117 2013-03-27 11:15:36Z suehring
52! Bugfix in OpenMP parallelization.
53!
54! 1113 2013-03-10 02:48:14Z raasch
55! GPU-porting of several loops, some loops rearranged
56!
57! 1111 2013-03-08 23:54:10Z
58! openACC statements added,
59! ibc_p_b = 2 removed
60!
61! 1092 2013-02-02 11:24:22Z raasch
62! unused variables removed
63!
64! 1036 2012-10-22 13:43:42Z raasch
65! code put under GPL (PALM 3.9)
66!
67! 1003 2012-09-14 14:35:53Z raasch
68! adjustment of array tend for cases with unequal subdomain sizes removed
69!
70! Revision 1.1  1997/07/24 11:24:44  raasch
71! Initial revision
72!
73!
74! Description:
75! ------------
76! Compute the divergence of the provisional velocity field. Solve the Poisson
77! equation for the perturbation pressure. Compute the final velocities using
78! this perturbation pressure. Compute the remaining divergence.
79!------------------------------------------------------------------------------!
80
81    USE arrays_3d,                                                             &
82        ONLY:  d, ddzu, ddzu_pres, ddzw, dzw, p, p_loc, tend, u, v, w
83
84    USE control_parameters,                                                    &
85        ONLY:  bc_lr_cyc, bc_ns_cyc, conserve_volume_flow, dt_3d,              &
86               gathered_size, ibc_p_b, ibc_p_t, intermediate_timestep_count,   &
87               mg_switch_to_pe0_level, on_device, outflow_l, outflow_n,        &
88               outflow_r, outflow_s, psolver, simulated_time, subdomain_size,  &
89               topography, volume_flow, volume_flow_area, volume_flow_initial
90
91    USE cpulog,                                                                &
92        ONLY:  cpu_log, log_point, log_point_s
93
94    USE grid_variables,                                                        &
95        ONLY:  ddx, ddy
96
97    USE indices,                                                               &
98        ONLY:  nbgp, ngp_2dh_outer, nx, nxl, nxlg, nxl_mg, nxr, nxrg, nxr_mg,  &
99               ny, nys, nysg, nys_mg, nyn, nyng, nyn_mg, nzb, nzb_s_inner,     &
100               nzb_u_inner, nzb_v_inner, nzb_w_inner, nzb_2d, nzt, nzt_mg,     &
101               rflags_s_inner
102
103    USE kinds
104
105    USE pegrid
106
107    USE poisfft_mod,                                                           &
108        ONLY:  poisfft
109
110    USE statistics,                                                            &
111        ONLY:  statistic_regions, sums_divnew_l, sums_divold_l, weight_pres,   &
112               weight_substep
113
114    IMPLICIT NONE
115
116    INTEGER(iwp) ::  i              !:
117    INTEGER(iwp) ::  j              !:
118    INTEGER(iwp) ::  k              !:
119
120    REAL(wp)     ::  ddt_3d         !:
121    REAL(wp)     ::  localsum       !:
122    REAL(wp)     ::  threadsum      !:
123    REAL(wp)     ::  d_weight_pres  !:
124
125    REAL(wp), DIMENSION(1:2)   ::  volume_flow_l       !:
126    REAL(wp), DIMENSION(1:2)   ::  volume_flow_offset  !:
127    REAL(wp), DIMENSION(1:nzt) ::  w_l                 !:
128    REAL(wp), DIMENSION(1:nzt) ::  w_l_l               !:
129
130
131    CALL cpu_log( log_point(8), 'pres', 'start' )
132
133
134    ddt_3d = 1.0 / dt_3d
135    d_weight_pres = 1.0 / weight_pres(intermediate_timestep_count)
136
137!
138!-- Multigrid method expects array d to have one ghost layer.
139!--
140    IF ( psolver == 'multigrid' )  THEN
141     
142       DEALLOCATE( d )
143       ALLOCATE( d(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ) 
144
145!
146!--    Since p is later used to hold the weighted average of the substeps, it
147!--    cannot be used in the iterative solver. Therefore, its initial value is
148!--    stored on p_loc, which is then iteratively advanced in every substep.
149       IF ( intermediate_timestep_count == 1 )  THEN
150          DO  i = nxl-1, nxr+1
151             DO  j = nys-1, nyn+1
152                DO  k = nzb, nzt+1
153                   p_loc(k,j,i) = p(k,j,i)
154                ENDDO
155             ENDDO
156          ENDDO
157       ENDIF
158       
159    ELSEIF ( psolver == 'sor'  .AND.  intermediate_timestep_count == 1 )  THEN
160
161!
162!--    Since p is later used to hold the weighted average of the substeps, it
163!--    cannot be used in the iterative solver. Therefore, its initial value is
164!--    stored on p_loc, which is then iteratively advanced in every substep.
165       p_loc = p
166
167    ENDIF
168
169!
170!-- Conserve the volume flow at the outflow in case of non-cyclic lateral
171!-- boundary conditions
172!-- WARNING: so far, this conservation does not work at the left/south
173!--          boundary if the topography at the inflow differs from that at the
174!--          outflow! For this case, volume_flow_area needs adjustment!
175!
176!-- Left/right
177    IF ( conserve_volume_flow  .AND.  ( outflow_l .OR. outflow_r ) )  THEN
178
179       volume_flow(1)   = 0.0
180       volume_flow_l(1) = 0.0
181
182       IF ( outflow_l )  THEN
183          i = 0
184       ELSEIF ( outflow_r )  THEN
185          i = nx+1
186       ENDIF
187
188       DO  j = nys, nyn
189!
190!--       Sum up the volume flow through the south/north boundary
191          DO  k = nzb_2d(j,i)+1, nzt
192             volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k)
193          ENDDO
194       ENDDO
195
196#if defined( __parallel )   
197       IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
198       CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, &
199                           MPI_SUM, comm1dy, ierr )   
200#else
201       volume_flow = volume_flow_l 
202#endif
203       volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) ) &
204                               / volume_flow_area(1)
205
206       DO  j = nysg, nyng
207          DO  k = nzb_2d(j,i)+1, nzt
208             u(k,j,i) = u(k,j,i) + volume_flow_offset(1)
209          ENDDO
210       ENDDO
211
212    ENDIF
213
214!
215!-- South/north
216    IF ( conserve_volume_flow  .AND.  ( outflow_n .OR. outflow_s ) )  THEN
217
218       volume_flow(2)   = 0.0
219       volume_flow_l(2) = 0.0
220
221       IF ( outflow_s )  THEN
222          j = 0
223       ELSEIF ( outflow_n )  THEN
224          j = ny+1
225       ENDIF
226
227       DO  i = nxl, nxr
228!
229!--       Sum up the volume flow through the south/north boundary
230          DO  k = nzb_2d(j,i)+1, nzt
231             volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k)
232          ENDDO
233       ENDDO
234
235#if defined( __parallel )   
236       IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
237       CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL, &
238                           MPI_SUM, comm1dx, ierr )   
239#else
240       volume_flow = volume_flow_l 
241#endif
242       volume_flow_offset(2) = ( volume_flow_initial(2) - volume_flow(2) )    &
243                               / volume_flow_area(2)
244
245       DO  i = nxlg, nxrg
246          DO  k = nzb_v_inner(j,i)+1, nzt
247             v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
248          ENDDO
249       ENDDO
250
251    ENDIF
252
253!
254!-- Remove mean vertical velocity
255    IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
256       IF ( simulated_time > 0.0 )  THEN ! otherwise nzb_w_inner not yet known
257          w_l = 0.0;  w_l_l = 0.0
258          DO  i = nxl, nxr
259             DO  j = nys, nyn
260                DO  k = nzb_w_inner(j,i)+1, nzt
261                   w_l_l(k) = w_l_l(k) + w(k,j,i)
262                ENDDO
263             ENDDO
264          ENDDO
265#if defined( __parallel )   
266          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
267          CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, &
268                              comm2d, ierr )
269#else
270          w_l = w_l_l 
271#endif
272          DO  k = 1, nzt
273             w_l(k) = w_l(k) / ngp_2dh_outer(k,0)
274          ENDDO
275          DO  i = nxlg, nxrg
276             DO  j = nysg, nyng
277                DO  k = nzb_w_inner(j,i)+1, nzt
278                   w(k,j,i) = w(k,j,i) - w_l(k)
279                ENDDO
280             ENDDO
281          ENDDO
282       ENDIF
283    ENDIF
284
285!
286!-- Compute the divergence of the provisional velocity field.
287    CALL cpu_log( log_point_s(1), 'divergence', 'start' )
288
289    IF ( psolver == 'multigrid' )  THEN
290       !$OMP PARALLEL DO SCHEDULE( STATIC )
291       DO  i = nxl-1, nxr+1
292          DO  j = nys-1, nyn+1
293             DO  k = nzb, nzt+1
294                d(k,j,i) = 0.0
295             ENDDO
296          ENDDO
297       ENDDO
298    ELSE
299       !$OMP PARALLEL DO SCHEDULE( STATIC )
300       !$acc kernels present( d )
301       DO  i = nxl, nxr
302          DO  j = nys, nyn
303             DO  k = nzb+1, nzt
304                d(k,j,i) = 0.0
305             ENDDO
306          ENDDO
307       ENDDO
308       !$acc end kernels
309    ENDIF
310
311    localsum  = 0.0
312    threadsum = 0.0
313
314#if defined( __ibm )
315    !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
316    !$OMP DO SCHEDULE( STATIC )
317    DO  i = nxl, nxr
318       DO  j = nys, nyn
319          DO  k = nzb_s_inner(j,i)+1, nzt
320             d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
321                          ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
322                          ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d      &
323                                                                * d_weight_pres
324          ENDDO
325!
326!--       Compute possible PE-sum of divergences for flow_statistics
327          DO  k = nzb_s_inner(j,i)+1, nzt
328             threadsum = threadsum + ABS( d(k,j,i) )
329          ENDDO
330
331       ENDDO
332    ENDDO
333
334    localsum = localsum + threadsum * dt_3d * &
335                          weight_pres(intermediate_timestep_count)
336
337    !$OMP END PARALLEL
338#else
339
340    !$OMP PARALLEL PRIVATE (i,j,k)
341    !$OMP DO SCHEDULE( STATIC )
342    !$acc kernels present( d, ddzw, rflags_s_inner, u, v, w )
343    !$acc loop collapse( 3 )
344    DO  i = nxl, nxr
345       DO  j = nys, nyn
346          DO  k = 1, nzt
347             d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx +             &
348                        ( v(k,j+1,i) - v(k,j,i) ) * ddy +               &
349                        ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d  &
350                        * d_weight_pres * rflags_s_inner(k,j,i)
351          ENDDO
352       ENDDO
353    ENDDO
354    !$acc end kernels
355    !$OMP END PARALLEL
356
357!
358!-- Compute possible PE-sum of divergences for flow_statistics
359    !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
360    !$OMP DO SCHEDULE( STATIC )
361    !$acc parallel loop collapse(3) present( d ) reduction(+:threadsum)
362    DO  i = nxl, nxr
363       DO  j = nys, nyn
364          DO  k = nzb+1, nzt
365             threadsum = threadsum + ABS( d(k,j,i) )
366          ENDDO
367       ENDDO
368    ENDDO
369    !$acc end parallel loop
370    localsum = localsum + threadsum * dt_3d * &
371                          weight_pres(intermediate_timestep_count)
372    !$OMP END PARALLEL
373#endif
374
375!
376!-- For completeness, set the divergence sum of all statistic regions to those
377!-- of the total domain
378    sums_divold_l(0:statistic_regions) = localsum
379
380    CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
381
382!
383!-- Compute the pressure perturbation solving the Poisson equation
384    IF ( psolver(1:7) == 'poisfft' )  THEN
385
386!
387!--    Solve Poisson equation via FFT and solution of tridiagonal matrices
388       IF ( psolver == 'poisfft' )  THEN
389
390          CALL poisfft( d )
391
392       ENDIF
393
394!
395!--    Store computed perturbation pressure and set boundary condition in
396!--    z-direction
397       !$OMP PARALLEL DO
398       !$acc kernels present( d, tend )
399       DO  i = nxl, nxr
400          DO  j = nys, nyn
401             DO  k = nzb+1, nzt
402                tend(k,j,i) = d(k,j,i)
403             ENDDO
404          ENDDO
405       ENDDO
406       !$acc end kernels
407
408!
409!--    Bottom boundary:
410!--    This condition is only required for internal output. The pressure
411!--    gradient (dp(nzb+1)-dp(nzb))/dz is not used anywhere else.
412       IF ( ibc_p_b == 1 )  THEN
413!
414!--       Neumann (dp/dz = 0)
415          !$OMP PARALLEL DO
416          !$acc kernels present( nzb_s_inner, tend )
417          DO  i = nxlg, nxrg
418             DO  j = nysg, nyng
419                tend(nzb_s_inner(j,i),j,i) = tend(nzb_s_inner(j,i)+1,j,i)
420             ENDDO
421          ENDDO
422          !$acc end kernels
423
424       ELSE
425!
426!--       Dirichlet
427          !$OMP PARALLEL DO
428          !$acc kernels present( tend )
429          DO  i = nxlg, nxrg
430             DO  j = nysg, nyng
431                tend(nzb_s_inner(j,i),j,i) = 0.0
432             ENDDO
433          ENDDO
434          !$acc end kernels
435
436       ENDIF
437
438!
439!--    Top boundary
440       IF ( ibc_p_t == 1 )  THEN
441!
442!--       Neumann
443          !$OMP PARALLEL DO
444          !$acc kernels present( tend )
445          DO  i = nxlg, nxrg
446             DO  j = nysg, nyng
447                tend(nzt+1,j,i) = tend(nzt,j,i)
448             ENDDO
449          ENDDO
450          !$acc end kernels
451
452       ELSE
453!
454!--       Dirichlet
455          !$OMP PARALLEL DO
456          !$acc kernels present( tend )
457          DO  i = nxlg, nxrg
458             DO  j = nysg, nyng
459                tend(nzt+1,j,i) = 0.0
460             ENDDO
461          ENDDO
462          !$acc end kernels
463
464       ENDIF
465
466!
467!--    Exchange boundaries for p
468       IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
469          on_device = .TRUE.         ! to be removed after complete porting
470       ELSE                          ! of ghost point exchange
471          !$acc update host( tend )
472       ENDIF
473       CALL exchange_horiz( tend, nbgp )
474       IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
475          on_device = .FALSE.        ! to be removed after complete porting
476       ELSE                          ! of ghost point exchange
477          !$acc update device( tend )
478       ENDIF
479     
480    ELSEIF ( psolver == 'sor' )  THEN
481
482!
483!--    Solve Poisson equation for perturbation pressure using SOR-Red/Black
484!--    scheme
485       CALL sor( d, ddzu_pres, ddzw, p_loc )
486       tend = p_loc
487
488    ELSEIF ( psolver == 'multigrid' )  THEN
489
490!
491!--    Solve Poisson equation for perturbation pressure using Multigrid scheme,
492!--    array tend is used to store the residuals, logical exchange_mg is used
493!--    to discern data exchange in multigrid ( 1 ghostpoint ) and normal grid
494!--    ( nbgp ghost points ).
495
496!--    If the number of grid points of the gathered grid, which is collected
497!--    on PE0, is larger than the number of grid points of an PE, than array
498!--    tend will be enlarged.
499       IF ( gathered_size > subdomain_size )  THEN
500          DEALLOCATE( tend )
501          ALLOCATE( tend(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg(          &
502                    mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1,&
503                    nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg(                   &
504                    mg_switch_to_pe0_level)+1) )
505       ENDIF
506
507       CALL poismg( tend )
508
509       IF ( gathered_size > subdomain_size )  THEN
510          DEALLOCATE( tend )
511          ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
512       ENDIF
513
514!
515!--    Restore perturbation pressure on tend because this array is used
516!--    further below to correct the velocity fields
517       DO  i = nxl-1, nxr+1
518          DO  j = nys-1, nyn+1
519             DO  k = nzb, nzt+1
520                tend(k,j,i) = p_loc(k,j,i)
521             ENDDO
522          ENDDO
523       ENDDO
524
525    ENDIF
526
527!
528!-- Store perturbation pressure on array p, used for pressure data output.
529!-- Ghost layers are added in the output routines (except sor-method: see below)
530    IF ( intermediate_timestep_count == 1 )  THEN
531       !$OMP PARALLEL PRIVATE (i,j,k)
532       !$OMP DO
533       !$acc kernels present( p, tend, weight_substep )
534       !$acc loop independent
535       DO  i = nxl-1, nxr+1
536          !$acc loop independent
537          DO  j = nys-1, nyn+1
538             !$acc loop independent
539             DO  k = nzb, nzt+1
540                p(k,j,i) = tend(k,j,i) * &
541                           weight_substep(intermediate_timestep_count)
542             ENDDO
543          ENDDO
544       ENDDO
545       !$acc end kernels
546       !$OMP END PARALLEL
547
548    ELSE 
549       !$OMP PARALLEL PRIVATE (i,j,k)
550       !$OMP DO
551       !$acc kernels present( p, tend, weight_substep )
552       !$acc loop independent
553       DO  i = nxl-1, nxr+1
554          !$acc loop independent
555          DO  j = nys-1, nyn+1
556             !$acc loop independent
557             DO  k = nzb, nzt+1
558                p(k,j,i) = p(k,j,i) + tend(k,j,i) * &
559                           weight_substep(intermediate_timestep_count)
560             ENDDO
561          ENDDO
562       ENDDO
563       !$acc end kernels
564       !$OMP END PARALLEL
565
566    ENDIF
567       
568!
569!-- SOR-method needs ghost layers for the next timestep
570    IF ( psolver == 'sor' )  CALL exchange_horiz( p, nbgp )
571
572!
573!-- Correction of the provisional velocities with the current perturbation
574!-- pressure just computed
575    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc .OR. bc_ns_cyc ) )  THEN
576       volume_flow_l(1) = 0.0
577       volume_flow_l(2) = 0.0
578    ENDIF
579
580    !$OMP PARALLEL PRIVATE (i,j,k)
581    !$OMP DO
582    !$acc kernels present( ddzu, nzb_u_inner, nzb_v_inner, nzb_w_inner, tend, u, v, w, weight_pres )
583    !$acc loop independent
584    DO  i = nxl, nxr   
585       !$acc loop independent
586       DO  j = nys, nyn
587          !$acc loop independent
588          DO  k = 1, nzt
589             IF ( k > nzb_w_inner(j,i) )  THEN
590                w(k,j,i) = w(k,j,i) - dt_3d *                                 &
591                           ( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1) *      &
592                           weight_pres(intermediate_timestep_count)
593             ENDIF
594          ENDDO
595          !$acc loop independent
596          DO  k = 1, nzt
597             IF ( k > nzb_u_inner(j,i) )  THEN
598                u(k,j,i) = u(k,j,i) - dt_3d *                                 &
599                           ( tend(k,j,i) - tend(k,j,i-1) ) * ddx *            &
600                           weight_pres(intermediate_timestep_count)
601             ENDIF
602          ENDDO
603          !$acc loop independent
604          DO  k = 1, nzt
605             IF ( k > nzb_v_inner(j,i) )  THEN
606                v(k,j,i) = v(k,j,i) - dt_3d *                                 &
607                           ( tend(k,j,i) - tend(k,j-1,i) ) * ddy *            &
608                           weight_pres(intermediate_timestep_count)
609             ENDIF
610          ENDDO                                                         
611
612       ENDDO
613    ENDDO
614    !$acc end kernels
615    !$OMP END PARALLEL
616
617!
618!-- Sum up the volume flow through the right and north boundary
619    IF ( conserve_volume_flow  .AND.  bc_lr_cyc  .AND.  bc_ns_cyc  .AND.  &
620         nxr == nx )  THEN
621
622       !$OMP PARALLEL PRIVATE (j,k)
623       !$OMP DO
624       DO  j = nys, nyn
625          !$OMP CRITICAL
626          DO  k = nzb_2d(j,nx) + 1, nzt
627             volume_flow_l(1) = volume_flow_l(1) + u(k,j,nx) * dzw(k)
628          ENDDO
629          !$OMP END CRITICAL
630       ENDDO
631       !$OMP END PARALLEL
632
633    ENDIF
634
635    IF ( conserve_volume_flow  .AND.  bc_ns_cyc  .AND.  bc_lr_cyc  .AND.  &
636         nyn == ny )  THEN
637
638       !$OMP PARALLEL PRIVATE (i,k)
639       !$OMP DO
640       DO  i = nxl, nxr
641          !$OMP CRITICAL
642          DO  k = nzb_2d(ny,i) + 1, nzt
643             volume_flow_l(2) = volume_flow_l(2) + v(k,ny,i) * dzw(k)
644           ENDDO
645          !$OMP END CRITICAL
646       ENDDO
647       !$OMP END PARALLEL
648
649    ENDIF
650   
651!
652!-- Conserve the volume flow
653    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc  .AND.  bc_ns_cyc ) )  THEN
654
655#if defined( __parallel )   
656       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
657       CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 2, MPI_REAL, &
658                           MPI_SUM, comm2d, ierr ) 
659#else
660       volume_flow = volume_flow_l 
661#endif   
662
663       volume_flow_offset = ( volume_flow_initial - volume_flow ) / &
664                            volume_flow_area
665
666       !$OMP PARALLEL PRIVATE (i,j,k)
667       !$OMP DO
668       DO  i = nxl, nxr
669          DO  j = nys, nyn
670             DO  k = nzb_u_inner(j,i) + 1, nzt
671                u(k,j,i) = u(k,j,i) + volume_flow_offset(1)
672             ENDDO
673             DO k = nzb_v_inner(j,i) + 1, nzt
674                v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
675             ENDDO
676          ENDDO
677       ENDDO
678
679       !$OMP END PARALLEL
680
681    ENDIF
682
683!
684!-- Exchange of boundaries for the velocities
685    IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
686       on_device = .TRUE.         ! to be removed after complete porting
687    ELSE                          ! of ghost point exchange
688       !$acc update host( u, v, w )
689    ENDIF
690    CALL exchange_horiz( u, nbgp )
691    CALL exchange_horiz( v, nbgp )
692    CALL exchange_horiz( w, nbgp )
693    IF ( numprocs == 1 )  THEN    ! workaround for single-core GPU runs
694       on_device = .FALSE.        ! to be removed after complete porting
695    ELSE                          ! of ghost point exchange
696       !$acc update device( u, v, w )
697    ENDIF
698
699!
700!-- Compute the divergence of the corrected velocity field,
701!-- a possible PE-sum is computed in flow_statistics
702    CALL cpu_log( log_point_s(1), 'divergence', 'start' )
703    sums_divnew_l = 0.0
704
705!
706!-- d must be reset to zero because it can contain nonzero values below the
707!-- topography
708    IF ( topography /= 'flat' )  d = 0.0
709
710    localsum  = 0.0
711    threadsum = 0.0
712
713    !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
714    !$OMP DO SCHEDULE( STATIC )
715#if defined( __ibm )
716    DO  i = nxl, nxr
717       DO  j = nys, nyn
718          DO  k = nzb_s_inner(j,i)+1, nzt
719             d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
720                        ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
721                        ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
722          ENDDO
723          DO  k = nzb+1, nzt
724             threadsum = threadsum + ABS( d(k,j,i) )
725          ENDDO
726       ENDDO
727    ENDDO
728#else
729    !$acc kernels present( d, ddzw, rflags_s_inner, u, v, w )
730    !$acc loop collapse( 3 )
731    DO  i = nxl, nxr
732       DO  j = nys, nyn
733          DO  k = 1, nzt
734                d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * ddx +   &
735                             ( v(k,j+1,i) - v(k,j,i) ) * ddy +   &
736                             ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) &
737                           ) * rflags_s_inner(k,j,i)
738          ENDDO
739       ENDDO
740    ENDDO
741    !$acc end kernels
742!
743!-- Compute possible PE-sum of divergences for flow_statistics
744    !$acc parallel loop collapse(3) present( d ) reduction(+:threadsum)
745    DO  i = nxl, nxr
746       DO  j = nys, nyn
747          DO  k = nzb+1, nzt
748             threadsum = threadsum + ABS( d(k,j,i) )
749          ENDDO
750       ENDDO
751    ENDDO
752    !$acc end parallel loop
753#endif
754
755    localsum = localsum + threadsum
756    !$OMP END PARALLEL
757
758!
759!-- For completeness, set the divergence sum of all statistic regions to those
760!-- of the total domain
761    sums_divnew_l(0:statistic_regions) = localsum
762
763    CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
764
765    CALL cpu_log( log_point(8), 'pres', 'stop' )
766
767
768 END SUBROUTINE pres
Note: See TracBrowser for help on using the repository browser.