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

Last change on this file since 2233 was 2233, checked in by suehring, 4 years ago

last commit documented

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