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

Last change on this file since 2106 was 2101, checked in by suehring, 8 years ago

last commit documented

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