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

Last change on this file since 2749 was 2718, checked in by maronga, 7 years ago

deleting of deprecated files; headers updated where needed

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