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

Last change on this file since 3775 was 3655, checked in by knoop, 6 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

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