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

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

Anelastic approximation implemented

  • Property svn:keywords set to Id
File size: 27.9 KB
RevLine 
[1682]1!> @file pres.f90
[2000]2!------------------------------------------------------------------------------!
[1036]3! This file is part of PALM.
4!
[2000]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.
[1036]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!
[1818]17! Copyright 1997-2016 Leibniz Universitaet Hannover
[2000]18!------------------------------------------------------------------------------!
[1036]19!
[484]20! Current revisions:
[1212]21! ------------------
[2037]22! Anelastic approximation implemented
[1932]23!
[1919]24! Former revisions:
25! -----------------
26! $Id: pres.f90 2037 2016-10-26 11:15:40Z knoop $
27!
[2001]28! 2000 2016-08-20 18:09:15Z knoop
29! Forced header and separation lines into 80 columns
30!
[1933]31! 1932 2016-06-10 12:09:21Z suehring
32! Initial version of purely vertical nesting introduced.
33!
[1932]34! 1931 2016-06-10 12:06:59Z suehring
35! Rename multigrid into multigrid_noopt and multigrid_fast into multigrid
36!
[1930]37! 1929 2016-06-09 16:25:25Z suehring
38! Bugfix: weight_substep for initial call, replace by local variable
39!
[1919]40! 1918 2016-05-27 14:35:57Z raasch
[1918]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
[1321]48!
[1909]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!
[1846]53! 1845 2016-04-08 08:29:13Z raasch
54! nzb_2d replace by nzb_u|v_inner
55!
[1800]56! 1799 2016-04-05 08:35:55Z gronemeier
57! Bugfix: excluded third dimension from horizontal volume flow calculation
58!
[1763]59! 1762 2016-02-25 12:31:13Z hellstea
60! Introduction of nested domain feature
61!
[1683]62! 1682 2015-10-07 23:56:08Z knoop
63! Code annotations made doxygen readable
64!
[1576]65! 1575 2015-03-27 09:56:27Z raasch
66! poismg_fast + respective module added, adjustments for psolver-queries
67!
[1343]68! 1342 2014-03-26 17:04:47Z kanani
69! REAL constants defined as wp-kind
70!
[1321]71! 1320 2014-03-20 08:40:49Z raasch
[1320]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
[708]79!
[1319]80! 1318 2014-03-17 13:35:16Z raasch
81! module interfaces removed
82!
[1307]83! 1306 2014-03-13 14:30:59Z raasch
84! second argument removed from routine poisfft
85!
[1258]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!
[1222]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!
[1213]94! 1212 2013-08-15 08:46:27Z raasch
95! call of poisfft_hybrid removed
96!
[1118]97! 1117 2013-03-27 11:15:36Z suehring
98! Bugfix in OpenMP parallelization.
99!
[1114]100! 1113 2013-03-10 02:48:14Z raasch
101! GPU-porting of several loops, some loops rearranged
102!
[1112]103! 1111 2013-03-08 23:54:10Z
104! openACC statements added,
105! ibc_p_b = 2 removed
106!
[1093]107! 1092 2013-02-02 11:24:22Z raasch
108! unused variables removed
109!
[1037]110! 1036 2012-10-22 13:43:42Z raasch
111! code put under GPL (PALM 3.9)
112!
[1004]113! 1003 2012-09-14 14:35:53Z raasch
114! adjustment of array tend for cases with unequal subdomain sizes removed
115!
[1]116! Revision 1.1  1997/07/24 11:24:44  raasch
117! Initial revision
118!
119!
120! Description:
121! ------------
[1682]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.
[1]125!------------------------------------------------------------------------------!
[1682]126 SUBROUTINE pres
127 
[1]128
[1320]129    USE arrays_3d,                                                             &
[2037]130        ONLY:  d, ddzu, ddzu_pres, ddzw, dzw, p, p_loc, rho_air, rho_air_zw,   &
131               tend, u, v, w
[1320]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,   &
[1908]136               intermediate_timestep_count_max, mg_switch_to_pe0_level,        &
[1933]137               nest_domain, on_device, outflow_l, outflow_n, outflow_r,        &
[1918]138               outflow_s, psolver, subdomain_size, topography, volume_flow,    &
139               volume_flow_area, volume_flow_initial
[1320]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,     &
[1845]150               nzb_u_inner, nzb_v_inner, nzb_w_inner, nzt, nzt_mg,             &
[1320]151               rflags_s_inner
152
153    USE kinds
154
[1]155    USE pegrid
[1933]156   
157    USE pmc_interface,                                                         &
158        ONLY:  nesting_mode 
[1]159
[1320]160    USE poisfft_mod,                                                           &
161        ONLY:  poisfft
162
[1575]163    USE poismg_mod
164
[1320]165    USE statistics,                                                            &
166        ONLY:  statistic_regions, sums_divnew_l, sums_divold_l, weight_pres,   &
167               weight_substep
168
[1]169    IMPLICIT NONE
170
[1682]171    INTEGER(iwp) ::  i              !<
172    INTEGER(iwp) ::  j              !<
173    INTEGER(iwp) ::  k              !<
[1]174
[1682]175    REAL(wp)     ::  ddt_3d         !<
[1918]176    REAL(wp)     ::  d_weight_pres  !<
[1682]177    REAL(wp)     ::  localsum       !<
178    REAL(wp)     ::  threadsum      !<
[1918]179    REAL(wp)     ::  weight_pres_l  !<
[1929]180    REAL(wp)     ::  weight_substep_l !<
[1]181
[1762]182    REAL(wp), DIMENSION(1:3)   ::  volume_flow_l       !<
183    REAL(wp), DIMENSION(1:3)   ::  volume_flow_offset  !<
[1682]184    REAL(wp), DIMENSION(1:nzt) ::  w_l                 !<
185    REAL(wp), DIMENSION(1:nzt) ::  w_l_l               !<
[1]186
[1933]187    LOGICAL :: nest_domain_nvn      !<
[1]188
[1933]189
[1]190    CALL cpu_log( log_point(8), 'pres', 'start' )
191
[85]192
[1918]193!
194!-- Calculate quantities to be used locally
[1342]195    ddt_3d = 1.0_wp / dt_3d
[1918]196    IF ( intermediate_timestep_count == 0 )  THEN
197!
198!--    If pres is called before initial time step
[1929]199       weight_pres_l    = 1.0_wp
200       d_weight_pres    = 1.0_wp
201       weight_substep_l = 1.0_wp
[1918]202    ELSE
[1929]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)
[1918]206    ENDIF
[85]207
[1]208!
[707]209!-- Multigrid method expects array d to have one ghost layer.
210!--
[1575]211    IF ( psolver(1:9) == 'multigrid' )  THEN
[667]212     
[1]213       DEALLOCATE( d )
[667]214       ALLOCATE( d(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ) 
[707]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.
[1762]220       IF ( intermediate_timestep_count <= 1 )  THEN
[707]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
[667]228       ENDIF
229       
[1918]230    ELSEIF ( psolver == 'sor'  .AND.  intermediate_timestep_count <= 1 )  THEN
[707]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
[1]238    ENDIF
239
240!
[75]241!-- Conserve the volume flow at the outflow in case of non-cyclic lateral
242!-- boundary conditions
[106]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
[709]248    IF ( conserve_volume_flow  .AND.  ( outflow_l .OR. outflow_r ) )  THEN
[680]249
[1342]250       volume_flow(1)   = 0.0_wp
251       volume_flow_l(1) = 0.0_wp
[106]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
[1845]262          DO  k = nzb_u_inner(j,i)+1, nzt
[667]263             volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k)
[106]264          ENDDO
265       ENDDO
266
267#if defined( __parallel )   
[680]268       IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
[106]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
[709]274       volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) ) &
[106]275                               / volume_flow_area(1)
276
[667]277       DO  j = nysg, nyng
[1845]278          DO  k = nzb_u_inner(j,i)+1, nzt
[106]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
[709]287    IF ( conserve_volume_flow  .AND.  ( outflow_n .OR. outflow_s ) )  THEN
[106]288
[1342]289       volume_flow(2)   = 0.0_wp
290       volume_flow_l(2) = 0.0_wp
[75]291
[106]292       IF ( outflow_s )  THEN
293          j = 0
294       ELSEIF ( outflow_n )  THEN
[75]295          j = ny+1
[106]296       ENDIF
297
298       DO  i = nxl, nxr
[75]299!
[106]300!--       Sum up the volume flow through the south/north boundary
[1845]301          DO  k = nzb_v_inner(j,i)+1, nzt
[667]302             volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k)
[75]303          ENDDO
[106]304       ENDDO
305
[75]306#if defined( __parallel )   
[680]307       IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
[75]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) )    &
[106]314                               / volume_flow_area(2)
[75]315
[667]316       DO  i = nxlg, nxrg
[709]317          DO  k = nzb_v_inner(j,i)+1, nzt
[106]318             v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
[75]319          ENDDO
[106]320       ENDDO
[75]321
322    ENDIF
323
[76]324!
[1762]325!-- Remove mean vertical velocity in case that Neumann conditions are
[1933]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.
[1918]330!-- This cannot be done before the first initial time step because ngp_2dh_outer
331!-- is not yet known then.
[1933]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 )        &
[1918]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)
[76]345             ENDDO
346          ENDDO
[1918]347       ENDDO
[76]348#if defined( __parallel )   
[1918]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 )
[76]352#else
[1918]353       w_l = w_l_l
[76]354#endif
[1918]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)
[76]362             ENDDO
363          ENDDO
[1918]364       ENDDO
[76]365    ENDIF
[75]366
367!
[1]368!-- Compute the divergence of the provisional velocity field.
369    CALL cpu_log( log_point_s(1), 'divergence', 'start' )
370
[1575]371    IF ( psolver(1:9) == 'multigrid' )  THEN
[1]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
[1342]376                d(k,j,i) = 0.0_wp
[1]377             ENDDO
378          ENDDO
379       ENDDO
380    ELSE
381       !$OMP PARALLEL DO SCHEDULE( STATIC )
[1111]382       !$acc kernels present( d )
[1003]383       DO  i = nxl, nxr
384          DO  j = nys, nyn
385             DO  k = nzb+1, nzt
[1342]386                d(k,j,i) = 0.0_wp
[1]387             ENDDO
388          ENDDO
389       ENDDO
[1111]390       !$acc end kernels
[1]391    ENDIF
392
[1342]393    localsum  = 0.0_wp
394    threadsum = 0.0_wp
[1]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
[2037]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 
[1]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
[1918]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
[1908]420    ENDIF
[1]421    !$OMP END PARALLEL
422#else
423
[1111]424    !$OMP PARALLEL PRIVATE (i,j,k)
425    !$OMP DO SCHEDULE( STATIC )
[1221]426    !$acc kernels present( d, ddzw, rflags_s_inner, u, v, w )
427    !$acc loop collapse( 3 )
[1111]428    DO  i = nxl, nxr
429       DO  j = nys, nyn
430          DO  k = 1, nzt
[2037]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)
[1]436          ENDDO
437       ENDDO
[1111]438    ENDDO
439    !$acc end kernels
440    !$OMP END PARALLEL
[1]441
442!
[1908]443!-- Compute possible PE-sum of divergences for flow_statistics. Carry out
444!-- computation only at last Runge-Kutta substep.
[1918]445    IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.  &
446         intermediate_timestep_count == 0 )  THEN
[1908]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
[1]455          ENDDO
456       ENDDO
[1908]457       !$acc end parallel loop
[1918]458       localsum = localsum + threadsum * dt_3d * weight_pres_l
[1908]459       !$OMP END PARALLEL
460    ENDIF
[1]461#endif
462
463!
464!-- For completeness, set the divergence sum of all statistic regions to those
465!-- of the total domain
[1918]466    IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.  &
467         intermediate_timestep_count == 0 )  THEN
[1908]468       sums_divold_l(0:statistic_regions) = localsum
[1918]469    ENDIF
[1]470
471    CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
472
473!
474!-- Compute the pressure perturbation solving the Poisson equation
[1575]475    IF ( psolver == 'poisfft' )  THEN
[1]476
477!
478!--    Solve Poisson equation via FFT and solution of tridiagonal matrices
[1575]479       CALL poisfft( d )
[1212]480
[1]481!
482!--    Store computed perturbation pressure and set boundary condition in
483!--    z-direction
484       !$OMP PARALLEL DO
[1113]485       !$acc kernels present( d, tend )
[1]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
[1113]493       !$acc end kernels
[1]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
[1113]503          !$acc kernels present( nzb_s_inner, tend )
[667]504          DO  i = nxlg, nxrg
505             DO  j = nysg, nyng
[1]506                tend(nzb_s_inner(j,i),j,i) = tend(nzb_s_inner(j,i)+1,j,i)
507             ENDDO
508          ENDDO
[1113]509          !$acc end kernels
[1]510
511       ELSE
512!
513!--       Dirichlet
514          !$OMP PARALLEL DO
[1113]515          !$acc kernels present( tend )
[667]516          DO  i = nxlg, nxrg
517             DO  j = nysg, nyng
[1342]518                tend(nzb_s_inner(j,i),j,i) = 0.0_wp
[1]519             ENDDO
520          ENDDO
[1113]521          !$acc end kernels
[1]522
523       ENDIF
524
525!
526!--    Top boundary
527       IF ( ibc_p_t == 1 )  THEN
528!
529!--       Neumann
530          !$OMP PARALLEL DO
[1113]531          !$acc kernels present( tend )
[667]532          DO  i = nxlg, nxrg
533             DO  j = nysg, nyng
[1]534                tend(nzt+1,j,i) = tend(nzt,j,i)
535             ENDDO
536          ENDDO
[1113]537          !$acc end kernels
[1]538
539       ELSE
540!
541!--       Dirichlet
542          !$OMP PARALLEL DO
[1113]543          !$acc kernels present( tend )
[667]544          DO  i = nxlg, nxrg
545             DO  j = nysg, nyng
[1342]546                tend(nzt+1,j,i) = 0.0_wp
[1]547             ENDDO
548          ENDDO
[1113]549          !$acc end kernels
[1]550
551       ENDIF
552
553!
554!--    Exchange boundaries for p
[1113]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
[667]560       CALL exchange_horiz( tend, nbgp )
[1113]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
[1]566     
567    ELSEIF ( psolver == 'sor' )  THEN
568
569!
570!--    Solve Poisson equation for perturbation pressure using SOR-Red/Black
571!--    scheme
[707]572       CALL sor( d, ddzu_pres, ddzw, p_loc )
573       tend = p_loc
[1]574
[1575]575    ELSEIF ( psolver(1:9) == 'multigrid' )  THEN
[1]576
577!
578!--    Solve Poisson equation for perturbation pressure using Multigrid scheme,
[667]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 ).
[778]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
[1575]594       IF ( psolver == 'multigrid' )  THEN
595          CALL poismg( tend )
596       ELSE
[1931]597          CALL poismg_noopt( tend )
[1575]598       ENDIF
[707]599
[778]600       IF ( gathered_size > subdomain_size )  THEN
601          DEALLOCATE( tend )
602          ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
603       ENDIF
604
[1]605!
606!--    Restore perturbation pressure on tend because this array is used
607!--    further below to correct the velocity fields
[707]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
[667]615
[1]616    ENDIF
617
618!
[707]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)
[1918]621    IF ( intermediate_timestep_count <= 1 )  THEN
[707]622       !$OMP PARALLEL PRIVATE (i,j,k)
623       !$OMP DO
[1929]624       !$acc kernels present( p, tend, weight_substep_l )
[1257]625       !$acc loop independent
[707]626       DO  i = nxl-1, nxr+1
[1257]627          !$acc loop independent
[707]628          DO  j = nys-1, nyn+1
[1257]629             !$acc loop independent
[707]630             DO  k = nzb, nzt+1
631                p(k,j,i) = tend(k,j,i) * &
[1929]632                           weight_substep_l
[673]633             ENDDO
[1]634          ENDDO
[707]635       ENDDO
[1113]636       !$acc end kernels
[707]637       !$OMP END PARALLEL
638
[1762]639    ELSEIF ( intermediate_timestep_count > 1 )  THEN
[707]640       !$OMP PARALLEL PRIVATE (i,j,k)
641       !$OMP DO
[1929]642       !$acc kernels present( p, tend, weight_substep_l )
[1257]643       !$acc loop independent
[707]644       DO  i = nxl-1, nxr+1
[1257]645          !$acc loop independent
[707]646          DO  j = nys-1, nyn+1
[1257]647             !$acc loop independent
[707]648             DO  k = nzb, nzt+1
649                p(k,j,i) = p(k,j,i) + tend(k,j,i) * &
[1929]650                           weight_substep_l
[673]651             ENDDO
652          ENDDO
[707]653       ENDDO
[1113]654       !$acc end kernels
[707]655       !$OMP END PARALLEL
656
657    ENDIF
[673]658       
[707]659!
660!-- SOR-method needs ghost layers for the next timestep
661    IF ( psolver == 'sor' )  CALL exchange_horiz( p, nbgp )
[682]662
[1]663!
664!-- Correction of the provisional velocities with the current perturbation
665!-- pressure just computed
[709]666    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc .OR. bc_ns_cyc ) )  THEN
[1342]667       volume_flow_l(1) = 0.0_wp
668       volume_flow_l(2) = 0.0_wp
[1]669    ENDIF
[707]670
[1]671    !$OMP PARALLEL PRIVATE (i,j,k)
672    !$OMP DO
[1918]673    !$acc kernels present( ddzu, nzb_u_inner, nzb_v_inner, nzb_w_inner, tend, u, v, w )
[1257]674    !$acc loop independent
[673]675    DO  i = nxl, nxr   
[1257]676       !$acc loop independent
[1]677       DO  j = nys, nyn
[1257]678          !$acc loop independent
[1113]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) *      &
[1918]683                           weight_pres_l
[1113]684             ENDIF
[1]685          ENDDO
[1257]686          !$acc loop independent
[1113]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 *            &
[1918]691                           weight_pres_l
[1113]692             ENDIF
[1]693          ENDDO
[1257]694          !$acc loop independent
[1113]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 *            &
[1918]699                           weight_pres_l
[1113]700             ENDIF
[673]701          ENDDO                                                         
[1]702
703       ENDDO
704    ENDDO
[1113]705    !$acc end kernels
[1]706    !$OMP END PARALLEL
[1113]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
[1845]717          DO  k = nzb_u_inner(j,nx) + 1, nzt
[1113]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
[1845]733          DO  k = nzb_v_inner(ny,i) + 1, nzt
[1113]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
[673]741   
[1]742!
743!-- Conserve the volume flow
[707]744    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc  .AND.  bc_ns_cyc ) )  THEN
[1]745
746#if defined( __parallel )   
[622]747       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[1]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
[1799]754       volume_flow_offset(1:2) = ( volume_flow_initial(1:2) - volume_flow(1:2) ) / &
755                            volume_flow_area(1:2)
[1]756
757       !$OMP PARALLEL PRIVATE (i,j,k)
758       !$OMP DO
759       DO  i = nxl, nxr
760          DO  j = nys, nyn
[667]761             DO  k = nzb_u_inner(j,i) + 1, nzt
762                u(k,j,i) = u(k,j,i) + volume_flow_offset(1)
[719]763             ENDDO
[1845]764             DO  k = nzb_v_inner(j,i) + 1, nzt
[667]765                v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
766             ENDDO
[1]767          ENDDO
768       ENDDO
[667]769
[1]770       !$OMP END PARALLEL
771
772    ENDIF
773
774!
775!-- Exchange of boundaries for the velocities
[1113]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
[667]781    CALL exchange_horiz( u, nbgp )
782    CALL exchange_horiz( v, nbgp )
783    CALL exchange_horiz( w, nbgp )
[1113]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
[1]789
790!
791!-- Compute the divergence of the corrected velocity field,
[1908]792!-- A possible PE-sum is computed in flow_statistics. Carry out computation
793!-- only at last Runge-Kutta step.
[1918]794    IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.  &
795         intermediate_timestep_count == 0 )  THEN
[1908]796       CALL cpu_log( log_point_s(1), 'divergence', 'start' )
797       sums_divnew_l = 0.0_wp
[1]798
799!
[1908]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
[1]803
[1908]804       localsum  = 0.0_wp
805       threadsum = 0.0_wp
[1]806
[1908]807       !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
808       !$OMP DO SCHEDULE( STATIC )
[1]809#if defined( __ibm )
[1908]810       DO  i = nxl, nxr
811          DO  j = nys, nyn
812             DO  k = nzb_s_inner(j,i)+1, nzt
[2037]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)
[1908]817             ENDDO
818             DO  k = nzb+1, nzt
819                threadsum = threadsum + ABS( d(k,j,i) )
820             ENDDO
[1]821          ENDDO
822       ENDDO
823#else
[1908]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
[2037]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)
[1908]834             ENDDO
[1113]835          ENDDO
836       ENDDO
[1908]837       !$acc end kernels
[1113]838!
[1908]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
[1]846          ENDDO
847       ENDDO
[1908]848       !$acc end parallel loop
[1]849#endif
[667]850
[1908]851       localsum = localsum + threadsum
852       !$OMP END PARALLEL
[1]853
854!
[1908]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
[1]858
[1908]859       CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
[1]860
[1908]861    ENDIF
862
[1]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.