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

Last change on this file since 1483 was 1343, checked in by kanani, 10 years ago

last commit documented

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