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

Last change on this file since 1614 was 1576, checked in by raasch, 9 years ago

last commit documented

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