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

Last change on this file since 3849 was 3849, checked in by knoop, 2 years ago

Bugfix: added proper OpenACC support to disturb_field

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