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

Last change on this file since 1450 was 1343, checked in by kanani, 11 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 24.0 KB
RevLine 
[1]1 SUBROUTINE pres
2
[1036]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!
[1310]17! Copyright 1997-2014 Leibniz Universitaet Hannover
[1036]18!--------------------------------------------------------------------------------!
19!
[484]20! Current revisions:
[1212]21! ------------------
[1343]22!
23!
[1321]24! Former revisions:
25! -----------------
26! $Id: pres.f90 1343 2014-03-26 17:07:58Z heinze $
27!
[1343]28! 1342 2014-03-26 17:04:47Z kanani
29! REAL constants defined as wp-kind
30!
[1321]31! 1320 2014-03-20 08:40:49Z raasch
[1320]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
[708]39!
[1319]40! 1318 2014-03-17 13:35:16Z raasch
41! module interfaces removed
42!
[1307]43! 1306 2014-03-13 14:30:59Z raasch
44! second argument removed from routine poisfft
45!
[1258]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!
[1222]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!
[1213]54! 1212 2013-08-15 08:46:27Z raasch
55! call of poisfft_hybrid removed
56!
[1118]57! 1117 2013-03-27 11:15:36Z suehring
58! Bugfix in OpenMP parallelization.
59!
[1114]60! 1113 2013-03-10 02:48:14Z raasch
61! GPU-porting of several loops, some loops rearranged
62!
[1112]63! 1111 2013-03-08 23:54:10Z
64! openACC statements added,
65! ibc_p_b = 2 removed
66!
[1093]67! 1092 2013-02-02 11:24:22Z raasch
68! unused variables removed
69!
[1037]70! 1036 2012-10-22 13:43:42Z raasch
71! code put under GPL (PALM 3.9)
72!
[1004]73! 1003 2012-09-14 14:35:53Z raasch
74! adjustment of array tend for cases with unequal subdomain sizes removed
75!
[1]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
[1320]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
[1]111    USE pegrid
112
[1320]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
[1]120    IMPLICIT NONE
121
[1320]122    INTEGER(iwp) ::  i              !:
123    INTEGER(iwp) ::  j              !:
124    INTEGER(iwp) ::  k              !:
[1]125
[1320]126    REAL(wp)     ::  ddt_3d         !:
127    REAL(wp)     ::  localsum       !:
128    REAL(wp)     ::  threadsum      !:
129    REAL(wp)     ::  d_weight_pres  !:
[1]130
[1320]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               !:
[1]135
136
137    CALL cpu_log( log_point(8), 'pres', 'start' )
138
[85]139
[1342]140    ddt_3d = 1.0_wp / dt_3d
141    d_weight_pres = 1.0_wp / weight_pres(intermediate_timestep_count)
[85]142
[1]143!
[707]144!-- Multigrid method expects array d to have one ghost layer.
145!--
[1]146    IF ( psolver == 'multigrid' )  THEN
[667]147     
[1]148       DEALLOCATE( d )
[667]149       ALLOCATE( d(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ) 
[707]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
[667]163       ENDIF
164       
[707]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
[1]173    ENDIF
174
175!
[75]176!-- Conserve the volume flow at the outflow in case of non-cyclic lateral
177!-- boundary conditions
[106]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
[709]183    IF ( conserve_volume_flow  .AND.  ( outflow_l .OR. outflow_r ) )  THEN
[680]184
[1342]185       volume_flow(1)   = 0.0_wp
186       volume_flow_l(1) = 0.0_wp
[106]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
[709]197          DO  k = nzb_2d(j,i)+1, nzt
[667]198             volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k)
[106]199          ENDDO
200       ENDDO
201
202#if defined( __parallel )   
[680]203       IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
[106]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
[709]209       volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) ) &
[106]210                               / volume_flow_area(1)
211
[667]212       DO  j = nysg, nyng
[709]213          DO  k = nzb_2d(j,i)+1, nzt
[106]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
[709]222    IF ( conserve_volume_flow  .AND.  ( outflow_n .OR. outflow_s ) )  THEN
[106]223
[1342]224       volume_flow(2)   = 0.0_wp
225       volume_flow_l(2) = 0.0_wp
[75]226
[106]227       IF ( outflow_s )  THEN
228          j = 0
229       ELSEIF ( outflow_n )  THEN
[75]230          j = ny+1
[106]231       ENDIF
232
233       DO  i = nxl, nxr
[75]234!
[106]235!--       Sum up the volume flow through the south/north boundary
[709]236          DO  k = nzb_2d(j,i)+1, nzt
[667]237             volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k)
[75]238          ENDDO
[106]239       ENDDO
240
[75]241#if defined( __parallel )   
[680]242       IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
[75]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) )    &
[106]249                               / volume_flow_area(2)
[75]250
[667]251       DO  i = nxlg, nxrg
[709]252          DO  k = nzb_v_inner(j,i)+1, nzt
[106]253             v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
[75]254          ENDDO
[106]255       ENDDO
[75]256
257    ENDIF
258
[76]259!
260!-- Remove mean vertical velocity
261    IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
[1342]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
[76]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 )   
[622]272          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[709]273          CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, &
274                              comm2d, ierr )
[76]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
[667]281          DO  i = nxlg, nxrg
282             DO  j = nysg, nyng
[76]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
[75]290
291!
[1]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
[1342]300                d(k,j,i) = 0.0_wp
[1]301             ENDDO
302          ENDDO
303       ENDDO
304    ELSE
305       !$OMP PARALLEL DO SCHEDULE( STATIC )
[1111]306       !$acc kernels present( d )
[1003]307       DO  i = nxl, nxr
308          DO  j = nys, nyn
309             DO  k = nzb+1, nzt
[1342]310                d(k,j,i) = 0.0_wp
[1]311             ENDDO
312          ENDDO
313       ENDDO
[1111]314       !$acc end kernels
[1]315    ENDIF
316
[1342]317    localsum  = 0.0_wp
318    threadsum = 0.0_wp
[1]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
[85]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 + &
[673]328                          ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k) ) * ddt_3d      &
329                                                                * d_weight_pres 
[1]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
[707]340    localsum = localsum + threadsum * dt_3d * &
341                          weight_pres(intermediate_timestep_count)
[693]342
[1]343    !$OMP END PARALLEL
344#else
345
[1111]346    !$OMP PARALLEL PRIVATE (i,j,k)
347    !$OMP DO SCHEDULE( STATIC )
[1221]348    !$acc kernels present( d, ddzw, rflags_s_inner, u, v, w )
349    !$acc loop collapse( 3 )
[1111]350    DO  i = nxl, nxr
351       DO  j = nys, nyn
352          DO  k = 1, nzt
[1221]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)
[1]357          ENDDO
358       ENDDO
[1111]359    ENDDO
360    !$acc end kernels
361    !$OMP END PARALLEL
[1]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 )
[1221]367    !$acc parallel loop collapse(3) present( d ) reduction(+:threadsum)
[1]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
[1257]375    !$acc end parallel loop
[707]376    localsum = localsum + threadsum * dt_3d * &
377                          weight_pres(intermediate_timestep_count)
[1]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
[1212]395
[1306]396          CALL poisfft( d )
[1113]397
[1]398       ENDIF
399
400!
401!--    Store computed perturbation pressure and set boundary condition in
402!--    z-direction
403       !$OMP PARALLEL DO
[1113]404       !$acc kernels present( d, tend )
[1]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
[1113]412       !$acc end kernels
[1]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
[1113]422          !$acc kernels present( nzb_s_inner, tend )
[667]423          DO  i = nxlg, nxrg
424             DO  j = nysg, nyng
[1]425                tend(nzb_s_inner(j,i),j,i) = tend(nzb_s_inner(j,i)+1,j,i)
426             ENDDO
427          ENDDO
[1113]428          !$acc end kernels
[1]429
430       ELSE
431!
432!--       Dirichlet
433          !$OMP PARALLEL DO
[1113]434          !$acc kernels present( tend )
[667]435          DO  i = nxlg, nxrg
436             DO  j = nysg, nyng
[1342]437                tend(nzb_s_inner(j,i),j,i) = 0.0_wp
[1]438             ENDDO
439          ENDDO
[1113]440          !$acc end kernels
[1]441
442       ENDIF
443
444!
445!--    Top boundary
446       IF ( ibc_p_t == 1 )  THEN
447!
448!--       Neumann
449          !$OMP PARALLEL DO
[1113]450          !$acc kernels present( tend )
[667]451          DO  i = nxlg, nxrg
452             DO  j = nysg, nyng
[1]453                tend(nzt+1,j,i) = tend(nzt,j,i)
454             ENDDO
455          ENDDO
[1113]456          !$acc end kernels
[1]457
458       ELSE
459!
460!--       Dirichlet
461          !$OMP PARALLEL DO
[1113]462          !$acc kernels present( tend )
[667]463          DO  i = nxlg, nxrg
464             DO  j = nysg, nyng
[1342]465                tend(nzt+1,j,i) = 0.0_wp
[1]466             ENDDO
467          ENDDO
[1113]468          !$acc end kernels
[1]469
470       ENDIF
471
472!
473!--    Exchange boundaries for p
[1113]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
[667]479       CALL exchange_horiz( tend, nbgp )
[1113]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
[1]485     
486    ELSEIF ( psolver == 'sor' )  THEN
487
488!
489!--    Solve Poisson equation for perturbation pressure using SOR-Red/Black
490!--    scheme
[707]491       CALL sor( d, ddzu_pres, ddzw, p_loc )
492       tend = p_loc
[1]493
494    ELSEIF ( psolver == 'multigrid' )  THEN
495
496!
497!--    Solve Poisson equation for perturbation pressure using Multigrid scheme,
[667]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 ).
[778]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
[1]513       CALL poismg( tend )
[707]514
[778]515       IF ( gathered_size > subdomain_size )  THEN
516          DEALLOCATE( tend )
517          ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
518       ENDIF
519
[1]520!
521!--    Restore perturbation pressure on tend because this array is used
522!--    further below to correct the velocity fields
[707]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
[667]530
[1]531    ENDIF
532
533!
[707]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
[1113]539       !$acc kernels present( p, tend, weight_substep )
[1257]540       !$acc loop independent
[707]541       DO  i = nxl-1, nxr+1
[1257]542          !$acc loop independent
[707]543          DO  j = nys-1, nyn+1
[1257]544             !$acc loop independent
[707]545             DO  k = nzb, nzt+1
546                p(k,j,i) = tend(k,j,i) * &
547                           weight_substep(intermediate_timestep_count)
[673]548             ENDDO
[1]549          ENDDO
[707]550       ENDDO
[1113]551       !$acc end kernels
[707]552       !$OMP END PARALLEL
553
554    ELSE 
555       !$OMP PARALLEL PRIVATE (i,j,k)
556       !$OMP DO
[1113]557       !$acc kernels present( p, tend, weight_substep )
[1257]558       !$acc loop independent
[707]559       DO  i = nxl-1, nxr+1
[1257]560          !$acc loop independent
[707]561          DO  j = nys-1, nyn+1
[1257]562             !$acc loop independent
[707]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)
[673]566             ENDDO
567          ENDDO
[707]568       ENDDO
[1113]569       !$acc end kernels
[707]570       !$OMP END PARALLEL
571
572    ENDIF
[673]573       
[707]574!
575!-- SOR-method needs ghost layers for the next timestep
576    IF ( psolver == 'sor' )  CALL exchange_horiz( p, nbgp )
[682]577
[1]578!
579!-- Correction of the provisional velocities with the current perturbation
580!-- pressure just computed
[709]581    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc .OR. bc_ns_cyc ) )  THEN
[1342]582       volume_flow_l(1) = 0.0_wp
583       volume_flow_l(2) = 0.0_wp
[1]584    ENDIF
[707]585
[1]586    !$OMP PARALLEL PRIVATE (i,j,k)
587    !$OMP DO
[1113]588    !$acc kernels present( ddzu, nzb_u_inner, nzb_v_inner, nzb_w_inner, tend, u, v, w, weight_pres )
[1257]589    !$acc loop independent
[673]590    DO  i = nxl, nxr   
[1257]591       !$acc loop independent
[1]592       DO  j = nys, nyn
[1257]593          !$acc loop independent
[1113]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
[1]600          ENDDO
[1257]601          !$acc loop independent
[1113]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
[1]608          ENDDO
[1257]609          !$acc loop independent
[1113]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
[673]616          ENDDO                                                         
[1]617
618       ENDDO
619    ENDDO
[1113]620    !$acc end kernels
[1]621    !$OMP END PARALLEL
[1113]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
[673]656   
[1]657!
658!-- Conserve the volume flow
[707]659    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc  .AND.  bc_ns_cyc ) )  THEN
[1]660
661#if defined( __parallel )   
[622]662       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[1]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
[667]676             DO  k = nzb_u_inner(j,i) + 1, nzt
677                u(k,j,i) = u(k,j,i) + volume_flow_offset(1)
[719]678             ENDDO
679             DO k = nzb_v_inner(j,i) + 1, nzt
[667]680                v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
681             ENDDO
[1]682          ENDDO
683       ENDDO
[667]684
[1]685       !$OMP END PARALLEL
686
687    ENDIF
688
689!
690!-- Exchange of boundaries for the velocities
[1113]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
[667]696    CALL exchange_horiz( u, nbgp )
697    CALL exchange_horiz( v, nbgp )
698    CALL exchange_horiz( w, nbgp )
[1113]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
[1]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' )
[1342]709    sums_divnew_l = 0.0_wp
[1]710
711!
712!-- d must be reset to zero because it can contain nonzero values below the
713!-- topography
[1342]714    IF ( topography /= 'flat' )  d = 0.0_wp
[1]715
[1342]716    localsum  = 0.0_wp
717    threadsum = 0.0_wp
[1]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
[1221]735    !$acc kernels present( d, ddzw, rflags_s_inner, u, v, w )
736    !$acc loop collapse( 3 )
[1]737    DO  i = nxl, nxr
738       DO  j = nys, nyn
[1113]739          DO  k = 1, nzt
[1221]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)
[1113]744          ENDDO
745       ENDDO
746    ENDDO
747    !$acc end kernels
748!
749!-- Compute possible PE-sum of divergences for flow_statistics
[1221]750    !$acc parallel loop collapse(3) present( d ) reduction(+:threadsum)
[1113]751    DO  i = nxl, nxr
752       DO  j = nys, nyn
[1221]753          DO  k = nzb+1, nzt
[1]754             threadsum = threadsum + ABS( d(k,j,i) )
755          ENDDO
756       ENDDO
757    ENDDO
[1257]758    !$acc end parallel loop
[1]759#endif
[667]760
[1]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.