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

Last change on this file since 3295 was 3183, checked in by suehring, 6 years ago

last commit documented

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