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

Last change on this file since 3611 was 3347, checked in by suehring, 6 years ago

Offline nesting revised and separated from large_scale_forcing_mod; Time-dependent synthetic turbulence generator; bugfixes in USM/LSM radiation model and init_pegrid

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