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

Last change on this file since 1930 was 1930, checked in by suehring, 8 years ago

last commit documented

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