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

Last change on this file since 3083 was 3016, checked in by Giersch, 6 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
Line 
1!> @file pres.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2018 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: pres.f90 3016 2018-05-09 10:53:37Z gronemeier $
27! Dollar sign added before Id
28!
29! 2696 2017-12-14 17:12:51Z kanani
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
34! Corrected "Former revisions" section
35!
36! 2696 2017-12-14 17:12:51Z kanani
37! Change in file header (GPL part)
38! poismg_noopt modularized (MS)
39!
40! 2298 2017-06-29 09:28:18Z raasch
41! comment changed + some formatting
42!
43! 2233 2017-05-30 18:08:54Z suehring
44!
45! 2232 2017-05-30 17:47:52Z suehring
46! Adjustments to new topography and surface concept
47!
48! 2118 2017-01-17 16:38:49Z raasch
49! OpenACC directives and related code removed
50!
51! 2073 2016-11-30 14:34:05Z raasch
52! openmp bugfix for calculation of new divergence
53!
54! 2037 2016-10-26 11:15:40Z knoop
55! Anelastic approximation implemented
56!
57! 2000 2016-08-20 18:09:15Z knoop
58! Forced header and separation lines into 80 columns
59!
60! 1932 2016-06-10 12:09:21Z suehring
61! Initial version of purely vertical nesting introduced.
62!
63! 1931 2016-06-10 12:06:59Z suehring
64! Rename multigrid into multigrid_noopt and multigrid_fast into multigrid
65!
66! 1929 2016-06-09 16:25:25Z suehring
67! Bugfix: weight_substep for initial call, replace by local variable
68!
69! 1918 2016-05-27 14:35:57Z raasch
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
77!
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!
82! 1845 2016-04-08 08:29:13Z raasch
83! nzb_2d replace by nzb_u|v_inner
84!
85! 1799 2016-04-05 08:35:55Z gronemeier
86! Bugfix: excluded third dimension from horizontal volume flow calculation
87!
88! 1762 2016-02-25 12:31:13Z hellstea
89! Introduction of nested domain feature
90!
91! 1682 2015-10-07 23:56:08Z knoop
92! Code annotations made doxygen readable
93!
94! 1575 2015-03-27 09:56:27Z raasch
95! poismg_fast + respective module added, adjustments for psolver-queries
96!
97! 1342 2014-03-26 17:04:47Z kanani
98! REAL constants defined as wp-kind
99!
100! 1320 2014-03-20 08:40:49Z raasch
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
108!
109! 1318 2014-03-17 13:35:16Z raasch
110! module interfaces removed
111!
112! 1306 2014-03-13 14:30:59Z raasch
113! second argument removed from routine poisfft
114!
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!
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!
123! 1212 2013-08-15 08:46:27Z raasch
124! call of poisfft_hybrid removed
125!
126! 1117 2013-03-27 11:15:36Z suehring
127! Bugfix in OpenMP parallelization.
128!
129! 1113 2013-03-10 02:48:14Z raasch
130! GPU-porting of several loops, some loops rearranged
131!
132! 1111 2013-03-08 23:54:10Z
133! openACC statements added,
134! ibc_p_b = 2 removed
135!
136! 1092 2013-02-02 11:24:22Z raasch
137! unused variables removed
138!
139! 1036 2012-10-22 13:43:42Z raasch
140! code put under GPL (PALM 3.9)
141!
142! 1003 2012-09-14 14:35:53Z raasch
143! adjustment of array tend for cases with unequal subdomain sizes removed
144!
145! Revision 1.1  1997/07/24 11:24:44  raasch
146! Initial revision
147!
148!
149! Description:
150! ------------
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.
154!------------------------------------------------------------------------------!
155 SUBROUTINE pres
156 
157
158    USE arrays_3d,                                                             &
159        ONLY:  d, ddzu, ddzu_pres, ddzw, dzw, p, p_loc, rho_air, rho_air_zw,   &
160               tend, u, v, w
161
162    USE control_parameters,                                                    &
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
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,  &
178               ny, nys, nysg, nys_mg, nyn, nyng, nyn_mg, nzb, nzt, nzt_mg,     &
179               wall_flags_0
180
181    USE kinds
182
183    USE pegrid
184   
185    USE pmc_interface,                                                         &
186        ONLY:  nesting_mode 
187
188    USE poisfft_mod,                                                           &
189        ONLY:  poisfft
190
191    USE poismg_mod
192
193    USE poismg_noopt_mod
194
195    USE statistics,                                                            &
196        ONLY:  statistic_regions, sums_divnew_l, sums_divold_l, weight_pres,   &
197               weight_substep
198
199    USE surface_mod,                                                           &
200        ONLY :  bc_h
201
202    IMPLICIT NONE
203
204    INTEGER(iwp) ::  i              !<
205    INTEGER(iwp) ::  j              !<
206    INTEGER(iwp) ::  k              !<
207    INTEGER(iwp) ::  m              !<
208
209    REAL(wp)     ::  ddt_3d         !<
210    REAL(wp)     ::  d_weight_pres  !<
211    REAL(wp)     ::  localsum       !<
212    REAL(wp)     ::  threadsum      !<
213    REAL(wp)     ::  weight_pres_l  !<
214    REAL(wp)     ::  weight_substep_l !<
215
216    REAL(wp), DIMENSION(1:3)   ::  volume_flow_l       !<
217    REAL(wp), DIMENSION(1:3)   ::  volume_flow_offset  !<
218    REAL(wp), DIMENSION(1:nzt) ::  w_l                 !<
219    REAL(wp), DIMENSION(1:nzt) ::  w_l_l               !<
220
221    LOGICAL :: nest_domain_nvn      !<
222
223
224    CALL cpu_log( log_point(8), 'pres', 'start' )
225
226
227!
228!-- Calculate quantities to be used locally
229    ddt_3d = 1.0_wp / dt_3d
230    IF ( intermediate_timestep_count == 0 )  THEN
231!
232!--    If pres is called before initial time step
233       weight_pres_l    = 1.0_wp
234       d_weight_pres    = 1.0_wp
235       weight_substep_l = 1.0_wp
236    ELSE
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)
240    ENDIF
241
242!
243!-- Multigrid method expects array d to have one ghost layer.
244!--
245    IF ( psolver(1:9) == 'multigrid' )  THEN
246     
247       DEALLOCATE( d )
248       ALLOCATE( d(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ) 
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.
254       IF ( intermediate_timestep_count <= 1 )  THEN
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
262       ENDIF
263       
264    ELSEIF ( psolver == 'sor'  .AND.  intermediate_timestep_count <= 1 )  THEN
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
272    ENDIF
273
274!
275!-- Conserve the volume flow at the outflow in case of non-cyclic lateral
276!-- boundary conditions
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
282    IF ( conserve_volume_flow  .AND.  ( outflow_l .OR. outflow_r ) )  THEN
283
284       volume_flow(1)   = 0.0_wp
285       volume_flow_l(1) = 0.0_wp
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
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                                            )
301          ENDDO
302       ENDDO
303
304#if defined( __parallel )   
305       IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
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
311       volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) ) &
312                               / volume_flow_area(1)
313
314       DO  j = nysg, nyng
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                                            )
320          ENDDO
321       ENDDO
322
323    ENDIF
324
325!
326!-- South/north
327    IF ( conserve_volume_flow  .AND.  ( outflow_n .OR. outflow_s ) )  THEN
328
329       volume_flow(2)   = 0.0_wp
330       volume_flow_l(2) = 0.0_wp
331
332       IF ( outflow_s )  THEN
333          j = 0
334       ELSEIF ( outflow_n )  THEN
335          j = ny+1
336       ENDIF
337
338       DO  i = nxl, nxr
339!
340!--       Sum up the volume flow through the south/north boundary
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                                            )
346          ENDDO
347       ENDDO
348
349#if defined( __parallel )   
350       IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
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) )    &
357                               / volume_flow_area(2)
358
359       DO  i = nxlg, nxrg
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                                            )
365          ENDDO
366       ENDDO
367
368    ENDIF
369
370!
371!-- Remove mean vertical velocity in case that Neumann conditions are
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.
376!-- This cannot be done before the first initial time step because ngp_2dh_outer
377!-- is not yet known then.
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 )        &
385    THEN
386       w_l = 0.0_wp;  w_l_l = 0.0_wp
387       DO  i = nxl, nxr
388          DO  j = nys, nyn
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                                            )
394             ENDDO
395          ENDDO
396       ENDDO
397#if defined( __parallel )   
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 )
401#else
402       w_l = w_l_l
403#endif
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
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                                            )
414             ENDDO
415          ENDDO
416       ENDDO
417    ENDIF
418
419!
420!-- Compute the divergence of the provisional velocity field.
421    CALL cpu_log( log_point_s(1), 'divergence', 'start' )
422
423    IF ( psolver(1:9) == 'multigrid' )  THEN
424       !$OMP PARALLEL DO SCHEDULE( STATIC ) PRIVATE (i,j,k)
425       DO  i = nxl-1, nxr+1
426          DO  j = nys-1, nyn+1
427             DO  k = nzb, nzt+1
428                d(k,j,i) = 0.0_wp
429             ENDDO
430          ENDDO
431       ENDDO
432    ELSE
433       !$OMP PARALLEL DO SCHEDULE( STATIC ) PRIVATE (i,j,k)
434       DO  i = nxl, nxr
435          DO  j = nys, nyn
436             DO  k = nzb+1, nzt
437                d(k,j,i) = 0.0_wp
438             ENDDO
439          ENDDO
440       ENDDO
441    ENDIF
442
443    localsum  = 0.0_wp
444    threadsum = 0.0_wp
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
451          DO  k = nzb+1, nzt
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)           &
456                        ) * ddt_3d * d_weight_pres                             &
457                                   * MERGE( 1.0_wp, 0.0_wp,                    &
458                                            BTEST( wall_flags_0(k,j,i), 0 )    &
459                                          )
460          ENDDO
461!
462!--       Compute possible PE-sum of divergences for flow_statistics
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                                          )
468          ENDDO
469
470       ENDDO
471    ENDDO
472
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
476    ENDIF
477    !$OMP END PARALLEL
478#else
479
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
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)           &
489                        ) * ddt_3d * d_weight_pres                             &
490                                   * MERGE( 1.0_wp, 0.0_wp,                    &
491                                            BTEST( wall_flags_0(k,j,i), 0 )    &
492                                          )     
493          ENDDO
494       ENDDO
495    ENDDO
496    !$OMP END PARALLEL
497
498!
499!-- Compute possible PE-sum of divergences for flow_statistics. Carry out
500!-- computation only at last Runge-Kutta substep.
501    IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.  &
502         intermediate_timestep_count == 0 )  THEN
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
510          ENDDO
511       ENDDO
512       localsum = localsum + threadsum * dt_3d * weight_pres_l
513       !$OMP END PARALLEL
514    ENDIF
515#endif
516
517!
518!-- For completeness, set the divergence sum of all statistic regions to those
519!-- of the total domain
520    IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.  &
521         intermediate_timestep_count == 0 )  THEN
522       sums_divold_l(0:statistic_regions) = localsum
523    ENDIF
524
525    CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
526
527!
528!-- Compute the pressure perturbation solving the Poisson equation
529    IF ( psolver == 'poisfft' )  THEN
530
531!
532!--    Solve Poisson equation via FFT and solution of tridiagonal matrices
533       CALL poisfft( d )
534
535!
536!--    Store computed perturbation pressure and set boundary condition in
537!--    z-direction
538       !$OMP PARALLEL DO PRIVATE (i,j,k)
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!
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
558             i = bc_h(0)%i(m)
559             j = bc_h(0)%j(m)
560             k = bc_h(0)%k(m)
561             tend(k-1,j,i) = tend(k,j,i)
562          ENDDO
563!
564!--       Downward facing
565          !$OMP PARALLEL DO PRIVATE( i, j, k )
566          DO  m = 1, bc_h(1)%ns
567             i = bc_h(1)%i(m)
568             j = bc_h(1)%j(m)
569             k = bc_h(1)%k(m)
570             tend(k+1,j,i) = tend(k,j,i)
571          ENDDO
572
573       ELSE
574!
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
580             i = bc_h(0)%i(m)
581             j = bc_h(0)%j(m)
582             k = bc_h(0)%k(m)
583             tend(k-1,j,i) = 0.0_wp
584          ENDDO
585!
586!--       Downward facing
587          !$OMP PARALLEL DO PRIVATE( i, j, k )
588          DO  m = 1, bc_h(1)%ns
589             i = bc_h(1)%i(m)
590             j = bc_h(1)%j(m)
591             k = bc_h(1)%k(m)
592             tend(k+1,j,i) = 0.0_wp
593          ENDDO
594
595       ENDIF
596
597!
598!--    Top boundary
599       IF ( ibc_p_t == 1 )  THEN
600!
601!--       Neumann
602          !$OMP PARALLEL DO PRIVATE (i,j,k)
603          DO  i = nxlg, nxrg
604             DO  j = nysg, nyng
605                tend(nzt+1,j,i) = tend(nzt,j,i)
606             ENDDO
607          ENDDO
608
609       ELSE
610!
611!--       Dirichlet
612          !$OMP PARALLEL DO PRIVATE (i,j,k)
613          DO  i = nxlg, nxrg
614             DO  j = nysg, nyng
615                tend(nzt+1,j,i) = 0.0_wp
616             ENDDO
617          ENDDO
618
619       ENDIF
620
621!
622!--    Exchange boundaries for p
623       CALL exchange_horiz( tend, nbgp )
624     
625    ELSEIF ( psolver == 'sor' )  THEN
626
627!
628!--    Solve Poisson equation for perturbation pressure using SOR-Red/Black
629!--    scheme
630       CALL sor( d, ddzu_pres, ddzw, p_loc )
631       tend = p_loc
632
633    ELSEIF ( psolver(1:9) == 'multigrid' )  THEN
634
635!
636!--    Solve Poisson equation for perturbation pressure using Multigrid scheme,
637!--    array tend is used to store the residuals.
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
650       IF ( psolver == 'multigrid' )  THEN
651          CALL poismg( tend )
652       ELSE
653          CALL poismg_noopt( tend )
654       ENDIF
655
656       IF ( gathered_size > subdomain_size )  THEN
657          DEALLOCATE( tend )
658          ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
659       ENDIF
660
661!
662!--    Restore perturbation pressure on tend because this array is used
663!--    further below to correct the velocity fields
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
671
672    ENDIF
673
674!
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)
677    IF ( intermediate_timestep_count <= 1 )  THEN
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) * &
684                           weight_substep_l
685             ENDDO
686          ENDDO
687       ENDDO
688       !$OMP END PARALLEL
689
690    ELSEIF ( intermediate_timestep_count > 1 )  THEN
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) * &
697                           weight_substep_l
698             ENDDO
699          ENDDO
700       ENDDO
701       !$OMP END PARALLEL
702
703    ENDIF
704       
705!
706!-- SOR-method needs ghost layers for the next timestep
707    IF ( psolver == 'sor' )  CALL exchange_horiz( p, nbgp )
708
709!
710!-- Correction of the provisional velocities with the current perturbation
711!-- pressure just computed
712    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc .OR. bc_ns_cyc ) )  THEN
713       volume_flow_l(1) = 0.0_wp
714       volume_flow_l(2) = 0.0_wp
715    ENDIF
716
717    !$OMP PARALLEL PRIVATE (i,j,k)
718    !$OMP DO
719    DO  i = nxl, nxr   
720       DO  j = nys, nyn
721
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                                            )
729          ENDDO
730
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                                            )
738          ENDDO
739
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                                            )
747          ENDDO                                                         
748
749       ENDDO
750    ENDDO
751    !$OMP END PARALLEL
752
753!
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!
765!-- Sum up the volume flow through the right and north boundary
766    IF ( conserve_volume_flow  .AND.  bc_lr_cyc  .AND.  bc_ns_cyc  .AND.       &
767         nxr == nx )  THEN
768
769       !$OMP PARALLEL PRIVATE (j,k)
770       !$OMP DO
771       DO  j = nys, nyn
772          !$OMP CRITICAL
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                                            )
778          ENDDO
779          !$OMP END CRITICAL
780       ENDDO
781       !$OMP END PARALLEL
782
783    ENDIF
784
785    IF ( conserve_volume_flow  .AND.  bc_ns_cyc  .AND.  bc_lr_cyc  .AND.       &
786         nyn == ny )  THEN
787
788       !$OMP PARALLEL PRIVATE (i,k)
789       !$OMP DO
790       DO  i = nxl, nxr
791          !$OMP CRITICAL
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                                            )
797           ENDDO
798          !$OMP END CRITICAL
799       ENDDO
800       !$OMP END PARALLEL
801
802    ENDIF
803   
804!
805!-- Conserve the volume flow
806    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc  .AND.  bc_ns_cyc ) )  THEN
807
808#if defined( __parallel )   
809       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
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
816       volume_flow_offset(1:2) = ( volume_flow_initial(1:2) - volume_flow(1:2) ) / &
817                            volume_flow_area(1:2)
818
819       !$OMP PARALLEL PRIVATE (i,j,k)
820       !$OMP DO
821       DO  i = nxl, nxr
822          DO  j = nys, nyn
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                                            )
828             ENDDO
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                                            )
834             ENDDO
835          ENDDO
836       ENDDO
837
838       !$OMP END PARALLEL
839
840    ENDIF
841
842!
843!-- Exchange of boundaries for the velocities
844    CALL exchange_horiz( u, nbgp )
845    CALL exchange_horiz( v, nbgp )
846    CALL exchange_horiz( w, nbgp )
847
848!
849!-- Compute the divergence of the corrected velocity field,
850!-- A possible PE-sum is computed in flow_statistics. Carry out computation
851!-- only at last Runge-Kutta step.
852    IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.  &
853         intermediate_timestep_count == 0 )  THEN
854       CALL cpu_log( log_point_s(1), 'divergence', 'start' )
855       sums_divnew_l = 0.0_wp
856
857!
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
861
862       localsum  = 0.0_wp
863       threadsum = 0.0_wp
864
865       !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
866#if defined( __ibm )
867       !$OMP DO SCHEDULE( STATIC )
868       DO  i = nxl, nxr
869          DO  j = nys, nyn
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                                 )
878             ENDDO
879             DO  k = nzb+1, nzt
880                threadsum = threadsum + ABS( d(k,j,i) )
881             ENDDO
882          ENDDO
883       ENDDO
884#else
885       !$OMP DO SCHEDULE( STATIC )
886       DO  i = nxl, nxr
887          DO  j = nys, nyn
888             DO  k = nzb+1, nzt
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)        &
893                           ) * MERGE( 1.0_wp, 0.0_wp,                          &
894                                   BTEST( wall_flags_0(k,j,i), 0 )             &
895                                    )
896             ENDDO
897          ENDDO
898       ENDDO
899!
900!--    Compute possible PE-sum of divergences for flow_statistics
901       !$OMP DO SCHEDULE( STATIC )
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
907          ENDDO
908       ENDDO
909#endif
910
911       localsum = localsum + threadsum
912       !$OMP END PARALLEL
913
914!
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
918
919       CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
920
921    ENDIF
922
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.