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

Last change on this file since 2606 was 2298, checked in by raasch, 7 years ago

write_binary is of type LOGICAL now, MPI2-related code removed, obsolete variables removed, sendrecv_in_background related parts removed, missing variable descriptions added

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