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

Last change on this file since 2037 was 2037, checked in by knoop, 5 years ago

Anelastic approximation implemented

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