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

Last change on this file since 4129 was 4015, checked in by raasch, 5 years ago

all reals changed to double precision in order to work with 32-bit working precision, otherwise calculated time intervals would mostly give zero; variable child_domain_nvn eliminated

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