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

Last change on this file since 2716 was 2716, checked in by kanani, 4 years ago

Correction of "Former revisions" section

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