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

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

Adjustments according new topography and surface-modelling concept implemented

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