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

Last change on this file since 2065 was 2038, checked in by knoop, 8 years ago

last commit documented

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