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

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

Merge of branch palm4u into trunk

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