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

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