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

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

last commit documented

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