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

Last change on this file since 1960 was 1933, checked in by hellstea, 8 years ago

last commit documented

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