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

Last change on this file since 3741 was 3655, checked in by knoop, 6 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

  • Property svn:keywords set to Id
File size: 33.9 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 3655 2019-01-07 16:51:22Z hellstea $
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!$ACC DATA COPYOUT(d, tend) &
239!$ACC COPY(u, v, w, p) &
240!$ACC COPYIN(rho_air, rho_air_zw, ddzu, ddzw, wall_flags_0)
241
242!$ACC DATA &
243!$ACC COPYIN(bc_h(0:1)) &
244!$ACC COPYIN(bc_h(0)%i(1:bc_h(0)%ns)) &
245!$ACC COPYIN(bc_h(0)%j(1:bc_h(0)%ns)) &
246!$ACC COPYIN(bc_h(0)%k(1:bc_h(0)%ns)) &
247!$ACC COPYIN(bc_h(1)%i(1:bc_h(1)%ns)) &
248!$ACC COPYIN(bc_h(1)%j(1:bc_h(1)%ns)) &
249!$ACC COPYIN(bc_h(1)%k(1:bc_h(1)%ns))
250
251
252!
253!-- Calculate quantities to be used locally
254    ddt_3d = 1.0_wp / dt_3d
255    IF ( intermediate_timestep_count == 0 )  THEN
256!
257!--    If pres is called before initial time step
258       weight_pres_l    = 1.0_wp
259       d_weight_pres    = 1.0_wp
260       weight_substep_l = 1.0_wp
261    ELSE
262       weight_pres_l    = weight_pres(intermediate_timestep_count)
263       d_weight_pres    = 1.0_wp / weight_pres(intermediate_timestep_count)
264       weight_substep_l = weight_substep(intermediate_timestep_count)
265    ENDIF
266
267!
268!-- Multigrid method expects array d to have one ghost layer.
269!--
270    IF ( psolver(1:9) == 'multigrid' )  THEN
271     
272       DEALLOCATE( d )
273       ALLOCATE( d(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ) 
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       IF ( intermediate_timestep_count <= 1 )  THEN
280          DO  i = nxl-1, nxr+1
281             DO  j = nys-1, nyn+1
282                DO  k = nzb, nzt+1
283                   p_loc(k,j,i) = p(k,j,i)
284                ENDDO
285             ENDDO
286          ENDDO
287       ENDIF
288       
289    ELSEIF ( psolver == 'sor'  .AND.  intermediate_timestep_count <= 1 )  THEN
290
291!
292!--    Since p is later used to hold the weighted average of the substeps, it
293!--    cannot be used in the iterative solver. Therefore, its initial value is
294!--    stored on p_loc, which is then iteratively advanced in every substep.
295       p_loc = p
296
297    ENDIF
298
299!
300!-- Conserve the volume flow at the outflow in case of non-cyclic lateral
301!-- boundary conditions
302!-- WARNING: so far, this conservation does not work at the left/south
303!--          boundary if the topography at the inflow differs from that at the
304!--          outflow! For this case, volume_flow_area needs adjustment!
305!
306!-- Left/right
307    IF ( conserve_volume_flow  .AND.  ( bc_radiation_l .OR.                    &
308                                        bc_radiation_r ) )  THEN
309
310       volume_flow(1)   = 0.0_wp
311       volume_flow_l(1) = 0.0_wp
312
313       IF ( bc_radiation_l )  THEN
314          i = 0
315       ELSEIF ( bc_radiation_r )  THEN
316          i = nx+1
317       ENDIF
318
319       DO  j = nys, nyn
320!
321!--       Sum up the volume flow through the south/north boundary
322          DO  k = nzb+1, nzt
323             volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzw(k)           &
324                                     * MERGE( 1.0_wp, 0.0_wp,                  &
325                                              BTEST( wall_flags_0(k,j,i), 1 )  &
326                                            )
327          ENDDO
328       ENDDO
329
330#if defined( __parallel )   
331       IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
332       CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, &
333                           MPI_SUM, comm1dy, ierr )   
334#else
335       volume_flow = volume_flow_l 
336#endif
337       volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) ) &
338                               / volume_flow_area(1)
339
340       DO  j = nysg, nyng
341          DO  k = nzb+1, nzt
342             u(k,j,i) = u(k,j,i) + volume_flow_offset(1)                       &
343                                     * MERGE( 1.0_wp, 0.0_wp,                  &
344                                              BTEST( wall_flags_0(k,j,i), 1 )  &
345                                            )
346          ENDDO
347       ENDDO
348
349    ENDIF
350
351!
352!-- South/north
353    IF ( conserve_volume_flow  .AND.  ( bc_radiation_n .OR. bc_radiation_s ) )  THEN
354
355       volume_flow(2)   = 0.0_wp
356       volume_flow_l(2) = 0.0_wp
357
358       IF ( bc_radiation_s )  THEN
359          j = 0
360       ELSEIF ( bc_radiation_n )  THEN
361          j = ny+1
362       ENDIF
363
364       DO  i = nxl, nxr
365!
366!--       Sum up the volume flow through the south/north boundary
367          DO  k = nzb+1, nzt
368             volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzw(k)           &
369                                     * MERGE( 1.0_wp, 0.0_wp,                  &
370                                              BTEST( wall_flags_0(k,j,i), 2 )  &
371                                            )
372          ENDDO
373       ENDDO
374
375#if defined( __parallel )   
376       IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
377       CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL, &
378                           MPI_SUM, comm1dx, ierr )   
379#else
380       volume_flow = volume_flow_l 
381#endif
382       volume_flow_offset(2) = ( volume_flow_initial(2) - volume_flow(2) )    &
383                               / volume_flow_area(2)
384
385       DO  i = nxlg, nxrg
386          DO  k = nzb+1, nzt
387             v(k,j,i) = v(k,j,i) + volume_flow_offset(2)                       &
388                                     * MERGE( 1.0_wp, 0.0_wp,                  &
389                                              BTEST( wall_flags_0(k,j,i), 2 )  &
390                                            )
391          ENDDO
392       ENDDO
393
394    ENDIF
395
396!
397!-- Remove mean vertical velocity in case that Neumann conditions are
398!-- used both at bottom and top boundary, and if not a nested domain in a
399!-- normal nesting run. In case of vertical nesting, this must be done.
400!-- Therefore an auxiliary logical variable child_domain_nvn is used here, and
401!-- nvn stands for non-vertical nesting.
402!-- This cannot be done before the first initial time step because ngp_2dh_outer
403!-- is not yet known then.
404    child_domain_nvn = child_domain
405    IF ( child_domain .AND. nesting_mode == 'vertical' )  THEN
406       child_domain_nvn = .FALSE.
407    ENDIF
408
409    IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1  .AND.  .NOT. nesting_offline  .AND. &
410         .NOT. child_domain_nvn  .AND. intermediate_timestep_count /= 0 )       &
411    THEN
412       w_l = 0.0_wp;  w_l_l = 0.0_wp
413       DO  i = nxl, nxr
414          DO  j = nys, nyn
415             DO  k = nzb+1, nzt
416                w_l_l(k) = w_l_l(k) + w(k,j,i)                                 &
417                                     * MERGE( 1.0_wp, 0.0_wp,                  &
418                                              BTEST( wall_flags_0(k,j,i), 3 )  &
419                                            )
420             ENDDO
421          ENDDO
422       ENDDO
423#if defined( __parallel )   
424       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
425       CALL MPI_ALLREDUCE( w_l_l(1), w_l(1), nzt, MPI_REAL, MPI_SUM, &
426                           comm2d, ierr )
427#else
428       w_l = w_l_l
429#endif
430       DO  k = 1, nzt
431          w_l(k) = w_l(k) / ngp_2dh_outer(k,0)
432       ENDDO
433       DO  i = nxlg, nxrg
434          DO  j = nysg, nyng
435             DO  k = nzb+1, nzt
436                w(k,j,i) = w(k,j,i) - w_l(k)                                   &
437                                     * MERGE( 1.0_wp, 0.0_wp,                  &
438                                              BTEST( wall_flags_0(k,j,i), 3 )  &
439                                            )
440             ENDDO
441          ENDDO
442       ENDDO
443    ENDIF
444
445!
446!-- Compute the divergence of the provisional velocity field.
447    CALL cpu_log( log_point_s(1), 'divergence', 'start' )
448
449    IF ( psolver(1:9) == 'multigrid' )  THEN
450       !$OMP PARALLEL DO SCHEDULE( STATIC ) PRIVATE (i,j,k)
451       DO  i = nxl-1, nxr+1
452          DO  j = nys-1, nyn+1
453             DO  k = nzb, nzt+1
454                d(k,j,i) = 0.0_wp
455             ENDDO
456          ENDDO
457       ENDDO
458    ELSE
459       !$OMP PARALLEL DO SCHEDULE( STATIC ) PRIVATE (i,j,k)
460       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
461       !$ACC PRESENT(d)
462       DO  i = nxl, nxr
463          DO  j = nys, nyn
464             DO  k = nzb+1, nzt
465                d(k,j,i) = 0.0_wp
466             ENDDO
467          ENDDO
468       ENDDO
469    ENDIF
470
471    localsum  = 0.0_wp
472    threadsum = 0.0_wp
473
474#if defined( __ibm )
475    !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
476    !$OMP DO SCHEDULE( STATIC )
477    DO  i = nxl, nxr
478       DO  j = nys, nyn
479          DO  k = nzb+1, nzt
480             d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx +       &
481                          ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy +       &
482                          ( w(k,j,i)   * rho_air_zw(k) -                       &
483                            w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k)           &
484                        ) * ddt_3d * d_weight_pres                             &
485                                   * MERGE( 1.0_wp, 0.0_wp,                    &
486                                            BTEST( wall_flags_0(k,j,i), 0 )    &
487                                          )
488          ENDDO
489!
490!--       Compute possible PE-sum of divergences for flow_statistics
491          DO  k = nzb+1, nzt
492             threadsum = threadsum + ABS( d(k,j,i) )                           &
493                                   * MERGE( 1.0_wp, 0.0_wp,                    &
494                                            BTEST( wall_flags_0(k,j,i), 0 )    &
495                                          )
496          ENDDO
497
498       ENDDO
499    ENDDO
500
501    IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.  &
502         intermediate_timestep_count == 0 )  THEN
503       localsum = localsum + threadsum * dt_3d * weight_pres_l
504    ENDIF
505    !$OMP END PARALLEL
506#else
507
508    !$OMP PARALLEL PRIVATE (i,j,k)
509    !$OMP DO SCHEDULE( STATIC )
510    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
511    !$ACC PRESENT(u, v, w, rho_air, rho_air_zw, ddzw, wall_flags_0) &
512    !$ACC PRESENT(d)
513    DO  i = nxl, nxr
514       DO  j = nys, nyn
515          DO  k = 1, nzt
516             d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx +       &
517                          ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy +       &
518                          ( w(k,j,i)   * rho_air_zw(k) -                       &
519                            w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k)           &
520                        ) * ddt_3d * d_weight_pres                             &
521                                   * MERGE( 1.0_wp, 0.0_wp,                    &
522                                            BTEST( wall_flags_0(k,j,i), 0 )    &
523                                          )     
524          ENDDO
525       ENDDO
526    ENDDO
527    !$OMP END PARALLEL
528
529!
530!-- Compute possible PE-sum of divergences for flow_statistics. Carry out
531!-- computation only at last Runge-Kutta substep.
532    IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.  &
533         intermediate_timestep_count == 0 )  THEN
534       !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
535       !$OMP DO SCHEDULE( STATIC )
536       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
537       !$ACC REDUCTION(+:threadsum) COPY(threadsum) &
538       !$ACC PRESENT(d)
539       DO  i = nxl, nxr
540          DO  j = nys, nyn
541             DO  k = nzb+1, nzt
542                threadsum = threadsum + ABS( d(k,j,i) )
543             ENDDO
544          ENDDO
545       ENDDO
546       localsum = localsum + threadsum * dt_3d * weight_pres_l
547       !$OMP END PARALLEL
548    ENDIF
549#endif
550
551!
552!-- For completeness, set the divergence sum of all statistic regions to those
553!-- of the total domain
554    IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.  &
555         intermediate_timestep_count == 0 )  THEN
556       sums_divold_l(0:statistic_regions) = localsum
557    ENDIF
558
559    CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
560
561!
562!-- Compute the pressure perturbation solving the Poisson equation
563    IF ( psolver == 'poisfft' )  THEN
564
565!
566!--    Solve Poisson equation via FFT and solution of tridiagonal matrices
567       CALL poisfft( d )
568
569!
570!--    Store computed perturbation pressure and set boundary condition in
571!--    z-direction
572       !$OMP PARALLEL DO PRIVATE (i,j,k)
573       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
574       !$ACC PRESENT(d, tend)
575       DO  i = nxl, nxr
576          DO  j = nys, nyn
577             DO  k = nzb+1, nzt
578                tend(k,j,i) = d(k,j,i)
579             ENDDO
580          ENDDO
581       ENDDO
582
583!
584!--    Bottom boundary:
585!--    This condition is only required for internal output. The pressure
586!--    gradient (dp(nzb+1)-dp(nzb))/dz is not used anywhere else.
587       IF ( ibc_p_b == 1 )  THEN
588!
589!--       Neumann (dp/dz = 0). Using surfae data type, first for non-natural
590!--       surfaces, then for natural and urban surfaces
591!--       Upward facing
592          !$OMP PARALLEL DO PRIVATE( i, j, k )
593          !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
594          !$ACC PRESENT(bc_h, tend)
595          DO  m = 1, bc_h(0)%ns
596             i = bc_h(0)%i(m)
597             j = bc_h(0)%j(m)
598             k = bc_h(0)%k(m)
599             tend(k-1,j,i) = tend(k,j,i)
600          ENDDO
601!
602!--       Downward facing
603          !$OMP PARALLEL DO PRIVATE( i, j, k )
604          !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
605          !$ACC PRESENT(bc_h, tend)
606          DO  m = 1, bc_h(1)%ns
607             i = bc_h(1)%i(m)
608             j = bc_h(1)%j(m)
609             k = bc_h(1)%k(m)
610             tend(k+1,j,i) = tend(k,j,i)
611          ENDDO
612
613       ELSE
614!
615!--       Dirichlet. Using surface data type, first for non-natural
616!--       surfaces, then for natural and urban surfaces
617!--       Upward facing
618          !$OMP PARALLEL DO PRIVATE( i, j, k )
619          DO  m = 1, bc_h(0)%ns
620             i = bc_h(0)%i(m)
621             j = bc_h(0)%j(m)
622             k = bc_h(0)%k(m)
623             tend(k-1,j,i) = 0.0_wp
624          ENDDO
625!
626!--       Downward facing
627          !$OMP PARALLEL DO PRIVATE( i, j, k )
628          DO  m = 1, bc_h(1)%ns
629             i = bc_h(1)%i(m)
630             j = bc_h(1)%j(m)
631             k = bc_h(1)%k(m)
632             tend(k+1,j,i) = 0.0_wp
633          ENDDO
634
635       ENDIF
636
637!
638!--    Top boundary
639       IF ( ibc_p_t == 1 )  THEN
640!
641!--       Neumann
642          !$OMP PARALLEL DO PRIVATE (i,j,k)
643          DO  i = nxlg, nxrg
644             DO  j = nysg, nyng
645                tend(nzt+1,j,i) = tend(nzt,j,i)
646             ENDDO
647          ENDDO
648
649       ELSE
650!
651!--       Dirichlet
652          !$OMP PARALLEL DO PRIVATE (i,j,k)
653          !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
654          !$ACC PRESENT(tend)
655          DO  i = nxlg, nxrg
656             DO  j = nysg, nyng
657                tend(nzt+1,j,i) = 0.0_wp
658             ENDDO
659          ENDDO
660
661       ENDIF
662
663!
664!--    Exchange boundaries for p
665       CALL exchange_horiz( tend, nbgp )
666     
667    ELSEIF ( psolver == 'sor' )  THEN
668
669!
670!--    Solve Poisson equation for perturbation pressure using SOR-Red/Black
671!--    scheme
672       CALL sor( d, ddzu_pres, ddzw, p_loc )
673       tend = p_loc
674
675    ELSEIF ( psolver(1:9) == 'multigrid' )  THEN
676
677!
678!--    Solve Poisson equation for perturbation pressure using Multigrid scheme,
679!--    array tend is used to store the residuals.
680
681!--    If the number of grid points of the gathered grid, which is collected
682!--    on PE0, is larger than the number of grid points of an PE, than array
683!--    tend will be enlarged.
684       IF ( gathered_size > subdomain_size )  THEN
685          DEALLOCATE( tend )
686          ALLOCATE( tend(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg(          &
687                    mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1,&
688                    nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg(                   &
689                    mg_switch_to_pe0_level)+1) )
690       ENDIF
691
692       IF ( psolver == 'multigrid' )  THEN
693          CALL poismg( tend )
694       ELSE
695          CALL poismg_noopt( tend )
696       ENDIF
697
698       IF ( gathered_size > subdomain_size )  THEN
699          DEALLOCATE( tend )
700          ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
701       ENDIF
702
703!
704!--    Restore perturbation pressure on tend because this array is used
705!--    further below to correct the velocity fields
706       DO  i = nxl-1, nxr+1
707          DO  j = nys-1, nyn+1
708             DO  k = nzb, nzt+1
709                tend(k,j,i) = p_loc(k,j,i)
710             ENDDO
711          ENDDO
712       ENDDO
713
714    ENDIF
715
716!
717!-- Store perturbation pressure on array p, used for pressure data output.
718!-- Ghost layers are added in the output routines (except sor-method: see below)
719    IF ( intermediate_timestep_count <= 1 )  THEN
720       !$OMP PARALLEL PRIVATE (i,j,k)
721       !$OMP DO
722       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
723       !$ACC PRESENT(p, tend)
724       DO  i = nxl-1, nxr+1
725          DO  j = nys-1, nyn+1
726             DO  k = nzb, nzt+1
727                p(k,j,i) = tend(k,j,i) * &
728                           weight_substep_l
729             ENDDO
730          ENDDO
731       ENDDO
732       !$OMP END PARALLEL
733
734    ELSEIF ( intermediate_timestep_count > 1 )  THEN
735       !$OMP PARALLEL PRIVATE (i,j,k)
736       !$OMP DO
737       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
738       !$ACC PRESENT(p, tend)
739       DO  i = nxl-1, nxr+1
740          DO  j = nys-1, nyn+1
741             DO  k = nzb, nzt+1
742                p(k,j,i) = p(k,j,i) + tend(k,j,i) * &
743                           weight_substep_l
744             ENDDO
745          ENDDO
746       ENDDO
747       !$OMP END PARALLEL
748
749    ENDIF
750       
751!
752!-- SOR-method needs ghost layers for the next timestep
753    IF ( psolver == 'sor' )  CALL exchange_horiz( p, nbgp )
754
755!
756!-- Correction of the provisional velocities with the current perturbation
757!-- pressure just computed
758    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc .OR. bc_ns_cyc ) )  THEN
759       volume_flow_l(1) = 0.0_wp
760       volume_flow_l(2) = 0.0_wp
761    ENDIF
762!
763!-- Add pressure gradients to the velocity components. Note, the loops are
764!-- running over the entire model domain, even though, in case of non-cyclic
765!-- boundaries u- and v-component are not prognostic at i=0 and j=0,
766!-- respectiveley. However, in case of Dirichlet boundary conditions for the
767!-- velocities, zero-gradient conditions for the pressure are set, so that
768!-- no modification is imposed at the boundaries.
769    !$OMP PARALLEL PRIVATE (i,j,k)
770    !$OMP DO
771    !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k) &
772    !$ACC PRESENT(u, v, w, tend, ddzu, wall_flags_0)
773    DO  i = nxl, nxr   
774       DO  j = nys, nyn
775
776          DO  k = nzb+1, nzt
777             w(k,j,i) = w(k,j,i) - dt_3d *                                     &
778                           ( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1)         &
779                                     * weight_pres_l                           &
780                                     * MERGE( 1.0_wp, 0.0_wp,                  &
781                                              BTEST( wall_flags_0(k,j,i), 3 )  &
782                                            )
783          ENDDO
784
785          DO  k = nzb+1, nzt
786             u(k,j,i) = u(k,j,i) - dt_3d *                                     &
787                           ( tend(k,j,i) - tend(k,j,i-1) ) * ddx               &
788                                     * weight_pres_l                           &
789                                     * MERGE( 1.0_wp, 0.0_wp,                  &
790                                              BTEST( wall_flags_0(k,j,i), 1 )  &
791                                            )
792          ENDDO
793
794          DO  k = nzb+1, nzt
795             v(k,j,i) = v(k,j,i) - dt_3d *                                     &
796                           ( tend(k,j,i) - tend(k,j-1,i) ) * ddy               &
797                                     * weight_pres_l                           &
798                                     * MERGE( 1.0_wp, 0.0_wp,                  &
799                                              BTEST( wall_flags_0(k,j,i), 2 )  &
800                                            )
801          ENDDO                                                         
802
803       ENDDO
804    ENDDO
805    !$OMP END PARALLEL
806
807!
808!-- The vertical velocity is not set to zero at nzt + 1 for nested domains
809!-- Instead it is set to the values of nzt (see routine vnest_boundary_conds
810!-- or pmci_interp_tril_t) BEFORE calling the pressure solver. To avoid jumps
811!-- while plotting profiles w at the top has to be set to the values in the
812!-- height nzt after above modifications. Hint: w level nzt+1 does not impact
813!-- results.
814    IF ( child_domain  .OR.  coupling_mode == 'vnested_fine' ) THEN
815       w(nzt+1,:,:) = w(nzt,:,:)
816    ENDIF
817
818!
819!-- Sum up the volume flow through the right and north boundary
820    IF ( conserve_volume_flow  .AND.  bc_lr_cyc  .AND.  bc_ns_cyc  .AND.       &
821         nxr == nx )  THEN
822
823       !$OMP PARALLEL PRIVATE (j,k)
824       !$OMP DO
825       DO  j = nys, nyn
826          !$OMP CRITICAL
827          DO  k = nzb+1, nzt
828             volume_flow_l(1) = volume_flow_l(1) + u(k,j,nxr) * dzw(k)         &
829                                     * MERGE( 1.0_wp, 0.0_wp,                  &
830                                              BTEST( wall_flags_0(k,j,nxr), 1 )&
831                                            )
832          ENDDO
833          !$OMP END CRITICAL
834       ENDDO
835       !$OMP END PARALLEL
836
837    ENDIF
838
839    IF ( conserve_volume_flow  .AND.  bc_ns_cyc  .AND.  bc_lr_cyc  .AND.       &
840         nyn == ny )  THEN
841
842       !$OMP PARALLEL PRIVATE (i,k)
843       !$OMP DO
844       DO  i = nxl, nxr
845          !$OMP CRITICAL
846          DO  k = nzb+1, nzt
847             volume_flow_l(2) = volume_flow_l(2) + v(k,nyn,i) * dzw(k)         &
848                                     * MERGE( 1.0_wp, 0.0_wp,                  &
849                                              BTEST( wall_flags_0(k,nyn,i), 2 )&
850                                            )
851           ENDDO
852          !$OMP END CRITICAL
853       ENDDO
854       !$OMP END PARALLEL
855
856    ENDIF
857   
858!
859!-- Conserve the volume flow
860    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc  .AND.  bc_ns_cyc ) )  THEN
861
862#if defined( __parallel )   
863       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
864       CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 2, MPI_REAL, &
865                           MPI_SUM, comm2d, ierr ) 
866#else
867       volume_flow = volume_flow_l 
868#endif   
869
870       volume_flow_offset(1:2) = ( volume_flow_initial(1:2) - volume_flow(1:2) ) / &
871                            volume_flow_area(1:2)
872
873       !$OMP PARALLEL PRIVATE (i,j,k)
874       !$OMP DO
875       DO  i = nxl, nxr
876          DO  j = nys, nyn
877             DO  k = nzb+1, nzt
878                u(k,j,i) = u(k,j,i) + volume_flow_offset(1)                    &
879                                     * MERGE( 1.0_wp, 0.0_wp,                  &
880                                              BTEST( wall_flags_0(k,j,i), 1 )  &
881                                            )
882             ENDDO
883             DO  k = nzb+1, nzt
884                v(k,j,i) = v(k,j,i) + volume_flow_offset(2)                    &
885                                     * MERGE( 1.0_wp, 0.0_wp,                  &
886                                              BTEST( wall_flags_0(k,j,i), 2 )  &
887                                            )
888             ENDDO
889          ENDDO
890       ENDDO
891
892       !$OMP END PARALLEL
893
894    ENDIF
895
896!
897!-- Exchange of boundaries for the velocities
898    CALL exchange_horiz( u, nbgp )
899    CALL exchange_horiz( v, nbgp )
900    CALL exchange_horiz( w, nbgp )
901
902!
903!-- Compute the divergence of the corrected velocity field,
904!-- A possible PE-sum is computed in flow_statistics. Carry out computation
905!-- only at last Runge-Kutta step.
906    IF ( intermediate_timestep_count == intermediate_timestep_count_max  .OR.  &
907         intermediate_timestep_count == 0 )  THEN
908       CALL cpu_log( log_point_s(1), 'divergence', 'start' )
909       sums_divnew_l = 0.0_wp
910
911!
912!--    d must be reset to zero because it can contain nonzero values below the
913!--    topography
914       IF ( topography /= 'flat' )  d = 0.0_wp
915
916       localsum  = 0.0_wp
917       threadsum = 0.0_wp
918
919       !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
920#if defined( __ibm )
921       !$OMP DO SCHEDULE( STATIC )
922       DO  i = nxl, nxr
923          DO  j = nys, nyn
924             DO  k = nzb+1, nzt
925             d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx +       &
926                          ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy +       &
927                          ( w(k,j,i)   * rho_air_zw(k) -                       &
928                            w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k)           &
929                        ) * MERGE( 1.0_wp, 0.0_wp,                             &
930                                   BTEST( wall_flags_0(k,j,i), 0 )             &
931                                 )
932             ENDDO
933             DO  k = nzb+1, nzt
934                threadsum = threadsum + ABS( d(k,j,i) )
935             ENDDO
936          ENDDO
937       ENDDO
938#else
939       !$OMP DO SCHEDULE( STATIC )
940       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
941       !$ACC PRESENT(u, v, w, rho_air, rho_air_zw, ddzw, wall_flags_0) &
942       !$ACC PRESENT(d)
943       DO  i = nxl, nxr
944          DO  j = nys, nyn
945             DO  k = nzb+1, nzt
946                d(k,j,i) = ( ( u(k,j,i+1) - u(k,j,i) ) * rho_air(k) * ddx +    &
947                             ( v(k,j+1,i) - v(k,j,i) ) * rho_air(k) * ddy +    &
948                             ( w(k,j,i)   * rho_air_zw(k) -                    &
949                               w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k)        &
950                           ) * MERGE( 1.0_wp, 0.0_wp,                          &
951                                   BTEST( wall_flags_0(k,j,i), 0 )             &
952                                    )
953             ENDDO
954          ENDDO
955       ENDDO
956!
957!--    Compute possible PE-sum of divergences for flow_statistics
958       !$OMP DO SCHEDULE( STATIC )
959       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
960       !$ACC REDUCTION(+:threadsum) COPY(threadsum) &
961       !$ACC PRESENT(d)
962       DO  i = nxl, nxr
963          DO  j = nys, nyn
964             DO  k = nzb+1, nzt
965                threadsum = threadsum + ABS( d(k,j,i) )
966             ENDDO
967          ENDDO
968       ENDDO
969#endif
970
971       localsum = localsum + threadsum
972       !$OMP END PARALLEL
973
974!
975!--    For completeness, set the divergence sum of all statistic regions to those
976!--    of the total domain
977       sums_divnew_l(0:statistic_regions) = localsum
978
979       CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
980
981    ENDIF
982
983    CALL cpu_log( log_point(8), 'pres', 'stop' )
984
985!$ACC END DATA
986!$ACC END DATA
987
988 END SUBROUTINE pres
Note: See TracBrowser for help on using the repository browser.