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

Last change on this file since 3016 was 3016, checked in by Giersch, 7 years ago

Dollar sign added before Id; Revised structure of reading svf data according to PALM coding standard: svf_code_field/len and fsvf removed, error messages PA0493 and PA0494 added, allocation status of output arrays checked

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