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

Last change on this file since 75 was 75, checked in by raasch, 15 years ago

preliminary update for changes concerning non-cyclic boundary conditions

  • Property svn:keywords set to Id
File size: 16.4 KB
Line 
1 SUBROUTINE pres
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6! Volume flow control for non-cyclic boundary conditions added (currently only
7! for the north boundary!!), 2nd+3rd argument removed from exchange horiz
8!
9! Former revisions:
10! -----------------
11! $Id: pres.f90 75 2007-03-22 09:54:05Z raasch $
12! RCS Log replace by Id keyword, revision history cleaned up
13!
14! Revision 1.25  2006/04/26 13:26:12  raasch
15! OpenMP optimization (+localsum, threadsum)
16!
17! Revision 1.1  1997/07/24 11:24:44  raasch
18! Initial revision
19!
20!
21! Description:
22! ------------
23! Compute the divergence of the provisional velocity field. Solve the Poisson
24! equation for the perturbation pressure. Compute the final velocities using
25! this perturbation pressure. Compute the remaining divergence.
26!------------------------------------------------------------------------------!
27
28    USE arrays_3d
29    USE constants
30    USE control_parameters
31    USE cpulog
32    USE grid_variables
33    USE indices
34    USE interfaces
35    USE pegrid
36    USE poisfft_mod
37    USE poisfft_hybrid_mod
38    USE statistics
39
40    IMPLICIT NONE
41
42    INTEGER ::  i, j, k, sr
43
44    REAL    ::  localsum, threadsum
45
46    REAL, DIMENSION(1:2) ::  volume_flow_l, volume_flow_offset
47
48
49    CALL cpu_log( log_point(8), 'pres', 'start' )
50
51!
52!-- Multigrid method needs additional grid points for the divergence array
53    IF ( psolver == 'multigrid' )  THEN
54       DEALLOCATE( d )
55       ALLOCATE( d(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
56    ENDIF
57
58!
59!-- Conserve the volume flow at the outflow in case of non-cyclic lateral
60!-- boundary conditions
61    IF ( conserve_volume_flow  .AND.  bc_ns == 'radiation/dirichlet')  THEN
62
63       volume_flow(2)   = 0.0
64       volume_flow_l(2) = 0.0
65
66       IF ( nyn == ny )  THEN
67          j = ny+1
68          DO  i = nxl, nxr
69!
70!--          Sum up the volume flow through the north boundary
71             DO  k = nzb_2d(j,i) + 1, nzt
72                volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzu(k)
73             ENDDO
74          ENDDO
75       ENDIF
76#if defined( __parallel )   
77       CALL MPI_ALLREDUCE( volume_flow_l(2), volume_flow(2), 1, MPI_REAL, &
78                           MPI_SUM, comm1dx, ierr )   
79#else
80       volume_flow = volume_flow_l 
81#endif
82       volume_flow_offset(2) = ( volume_flow_initial(2) - volume_flow(2) )    &
83                               / volume_flow_area(2)                         
84
85       IF ( outflow_n )  THEN
86          j = nyn+1
87          DO  i = nxl, nxr
88             DO  k = nzb_v_inner(j,i) + 1, nzt
89                v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
90             ENDDO
91          ENDDO
92       ENDIF
93
94       CALL exchange_horiz( v )
95
96    ENDIF
97
98
99!
100!-- Compute the divergence of the provisional velocity field.
101    CALL cpu_log( log_point_s(1), 'divergence', 'start' )
102
103    IF ( psolver == 'multigrid' )  THEN
104       !$OMP PARALLEL DO SCHEDULE( STATIC )
105       DO  i = nxl-1, nxr+1
106          DO  j = nys-1, nyn+1
107             DO  k = nzb, nzt+1
108                d(k,j,i) = 0.0
109             ENDDO
110          ENDDO
111       ENDDO
112    ELSE
113       !$OMP PARALLEL DO SCHEDULE( STATIC )
114       DO  i = nxl, nxra
115          DO  j = nys, nyna
116             DO  k = nzb+1, nzta
117                d(k,j,i) = 0.0
118             ENDDO
119          ENDDO
120       ENDDO
121    ENDIF
122
123    localsum  = 0.0
124    threadsum = 0.0
125
126#if defined( __ibm )
127    !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
128    !$OMP DO SCHEDULE( STATIC )
129    DO  i = nxl, nxr
130       DO  j = nys, nyn
131          DO  k = nzb_s_inner(j,i)+1, nzt
132             d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
133                        ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
134                        ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
135          ENDDO
136!
137!--       Additional pressure boundary condition at the bottom boundary for
138!--       inhomogeneous Prandtl layer heat fluxes and temperatures, respectively
139!--       dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0.
140!--       This condition must not be applied at the start of a run, because then
141!--       flow_statistics has not yet been called and thus sums = 0.
142          IF ( ibc_p_b == 2  .AND.  sums(nzb+1,4) /= 0.0 )  THEN
143             k = nzb_s_inner(j,i)
144             d(k+1,j,i) = d(k+1,j,i) + (                                     &
145                                         ( usws(j,i+1) - usws(j,i) ) * ddx   &
146                                       + ( vsws(j+1,i) - vsws(j,i) ) * ddy   &
147                                       - g * ( pt(k+1,j,i) - sums(k+1,4) ) / &
148                                         sums(k+1,4)                         &
149                                       ) * ddzw(k+1)
150          ENDIF
151
152!
153!--       Compute possible PE-sum of divergences for flow_statistics
154          DO  k = nzb_s_inner(j,i)+1, nzt
155             threadsum = threadsum + ABS( d(k,j,i) )
156          ENDDO
157
158!
159!--       Velocity corrections are made with Euler step size. Right hand side
160!--       of Poisson equation has to be set appropriately
161          DO  k = nzb_s_inner(j,i)+1, nzt
162             d(k,j,i) = d(k,j,i) / dt_3d
163          ENDDO
164
165       ENDDO
166    ENDDO
167
168    localsum = localsum + threadsum
169    !$OMP END PARALLEL
170#else
171    IF ( ibc_p_b == 2 .AND. sums(nzb+1,4) /= 0.0 )  THEN
172       !$OMP PARALLEL PRIVATE (i,j,k)
173       !$OMP DO SCHEDULE( STATIC )
174       DO  i = nxl, nxr
175          DO  j = nys, nyn
176             DO  k = nzb_s_inner(j,i)+1, nzt
177                d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
178                           ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
179                           ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
180             ENDDO
181          ENDDO
182!
183!--       Additional pressure boundary condition at the bottom boundary for
184!--       inhomogeneous Prandtl layer heat fluxes and temperatures, respectively
185!--       dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0.
186!--       This condition must not be applied at the start of a run, because then
187!--       flow_statistics has not yet been called and thus sums = 0.
188          DO  j = nys, nyn
189              k = nzb_s_inner(j,i)
190              d(k+1,j,i) = d(k+1,j,i) + (                        &
191                             ( usws(j,i+1) - usws(j,i) ) * ddx   &
192                           + ( vsws(j+1,i) - vsws(j,i) ) * ddy   &
193                           - g * ( pt(k+1,j,i) - sums(k+1,4) ) / &
194                             sums(k+1,4)                         &
195                                        ) * ddzw(k+1)
196          ENDDO
197       ENDDO
198       !$OMP END PARALLEL
199
200    ELSE
201
202       !$OMP PARALLEL PRIVATE (i,j,k)
203       !$OMP DO SCHEDULE( STATIC )
204       DO  i = nxl, nxr
205          DO  j = nys, nyn
206             DO  k = nzb_s_inner(j,i)+1, nzt
207                d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
208                           ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
209                           ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
210             ENDDO
211          ENDDO
212       ENDDO
213       !$OMP END PARALLEL
214
215    ENDIF
216
217!
218!-- Compute possible PE-sum of divergences for flow_statistics
219    !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
220    !$OMP DO SCHEDULE( STATIC )
221    DO  i = nxl, nxr
222       DO  j = nys, nyn
223          DO  k = nzb+1, nzt
224             threadsum = threadsum + ABS( d(k,j,i) )
225          ENDDO
226       ENDDO
227    ENDDO
228    localsum = localsum + threadsum
229    !$OMP END PARALLEL
230
231!
232!-- Velocity corrections are made with Euler step size. Right hand side
233!-- of Poisson equation has to be set appropriately
234    !$OMP DO SCHEDULE( STATIC )
235    DO  i = nxl, nxr
236       DO  j = nys, nyn
237          DO  k = nzb_s_inner(j,i)+1, nzt
238             d(k,j,i) = d(k,j,i) / dt_3d
239          ENDDO
240       ENDDO
241    ENDDO
242#endif
243
244!
245!-- For completeness, set the divergence sum of all statistic regions to those
246!-- of the total domain
247    sums_divold_l(0:statistic_regions) = localsum
248
249!
250!-- Determine absolute minimum/maximum (only for test cases, therefore as
251!-- comment line)
252!    CALL global_min_max( nzb+1, nzt, nys, nyn, nxl, nxr, d, 'abs', divmax, &
253!                         divmax_ijk )
254
255    CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
256
257!
258!-- Compute the pressure perturbation solving the Poisson equation
259    IF ( psolver(1:7) == 'poisfft' )  THEN
260
261!
262!--    Enlarge the size of tend, used as a working array for the transpositions
263       IF ( nxra > nxr  .OR.  nyna > nyn  .OR.  nza > nz )  THEN
264          DEALLOCATE( tend )
265          ALLOCATE( tend(1:nza,nys:nyna,nxl:nxra) )
266       ENDIF
267
268!
269!--    Solve Poisson equation via FFT and solution of tridiagonal matrices
270       IF ( psolver == 'poisfft' )  THEN
271!
272!--       Solver for 2d-decomposition
273          CALL poisfft( d, tend )
274       ELSEIF ( psolver == 'poisfft_hybrid' )  THEN 
275!
276!--       Solver for 1d-decomposition (using MPI and OpenMP).
277!--       The old hybrid-solver is still included here, as long as there
278!--       are some optimization problems in poisfft
279          CALL poisfft_hybrid( d )
280       ENDIF
281
282!
283!--    Resize tend to its normal size
284       IF ( nxra > nxr  .OR.  nyna > nyn  .OR.  nza > nz )  THEN
285          DEALLOCATE( tend )
286          ALLOCATE( tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
287       ENDIF
288
289!
290!--    Store computed perturbation pressure and set boundary condition in
291!--    z-direction
292       !$OMP PARALLEL DO
293       DO  i = nxl, nxr
294          DO  j = nys, nyn
295             DO  k = nzb+1, nzt
296                tend(k,j,i) = d(k,j,i)
297             ENDDO
298          ENDDO
299       ENDDO
300
301!
302!--    Bottom boundary:
303!--    This condition is only required for internal output. The pressure
304!--    gradient (dp(nzb+1)-dp(nzb))/dz is not used anywhere else.
305       IF ( ibc_p_b == 1 )  THEN
306!
307!--       Neumann (dp/dz = 0)
308          !$OMP PARALLEL DO
309          DO  i = nxl-1, nxr+1
310             DO  j = nys-1, nyn+1
311                tend(nzb_s_inner(j,i),j,i) = tend(nzb_s_inner(j,i)+1,j,i)
312             ENDDO
313          ENDDO
314
315       ELSEIF ( ibc_p_b == 2 )  THEN
316!
317!--       Neumann condition for inhomogeneous surfaces,
318!--       here currently still in the form of a zero gradient. Actually
319!--       dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0 would have to be used for
320!--       the computation (cf. above: computation of divergences).
321          !$OMP PARALLEL DO
322          DO  i = nxl-1, nxr+1
323             DO  j = nys-1, nyn+1
324                tend(nzb_s_inner(j,i),j,i) = tend(nzb_s_inner(j,i)+1,j,i)
325             ENDDO
326          ENDDO
327
328       ELSE
329!
330!--       Dirichlet
331          !$OMP PARALLEL DO
332          DO  i = nxl-1, nxr+1
333             DO  j = nys-1, nyn+1
334                tend(nzb_s_inner(j,i),j,i) = 0.0
335             ENDDO
336          ENDDO
337
338       ENDIF
339
340!
341!--    Top boundary
342       IF ( ibc_p_t == 1 )  THEN
343!
344!--       Neumann
345          !$OMP PARALLEL DO
346          DO  i = nxl-1, nxr+1
347             DO  j = nys-1, nyn+1
348                tend(nzt+1,j,i) = tend(nzt,j,i)
349             ENDDO
350          ENDDO
351
352       ELSE
353!
354!--       Dirichlet
355          !$OMP PARALLEL DO
356          DO  i = nxl-1, nxr+1
357             DO  j = nys-1, nyn+1
358                tend(nzt+1,j,i) = 0.0
359             ENDDO
360          ENDDO
361
362       ENDIF
363
364!
365!--    Exchange boundaries for p
366       CALL exchange_horiz( tend )
367     
368    ELSEIF ( psolver == 'sor' )  THEN
369
370!
371!--    Solve Poisson equation for perturbation pressure using SOR-Red/Black
372!--    scheme
373       CALL sor( d, ddzu, ddzw, p )
374       tend = p
375
376    ELSEIF ( psolver == 'multigrid' )  THEN
377
378!
379!--    Solve Poisson equation for perturbation pressure using Multigrid scheme,
380!--    array tend is used to store the residuals
381       CALL poismg( tend )
382 
383!
384!--    Restore perturbation pressure on tend because this array is used
385!--    further below to correct the velocity fields
386       tend = p
387
388    ENDIF
389
390!
391!-- Store perturbation pressure on array p, used in the momentum equations
392    IF ( psolver(1:7) == 'poisfft' )  THEN
393!
394!--    Here, only the values from the left and right boundaries are copied
395!--    The remaining values are copied in the following loop due to speed
396!--    optimization
397       !$OMP PARALLEL DO
398       DO  j = nys-1, nyn+1
399          DO  k = nzb, nzt+1
400             p(k,j,nxl-1) = tend(k,j,nxl-1)
401             p(k,j,nxr+1) = tend(k,j,nxr+1)
402          ENDDO
403       ENDDO
404    ENDIF
405
406!
407!-- Correction of the provisional velocities with the current perturbation
408!-- pressure just computed
409    IF ( conserve_volume_flow  .AND. &
410         ( bc_lr == 'cyclic'  .OR.  bc_ns == 'cyclic' ) )  THEN
411       volume_flow_l(1) = 0.0
412       volume_flow_l(2) = 0.0
413    ENDIF
414    !$OMP PARALLEL PRIVATE (i,j,k)
415    !$OMP DO
416    DO  i = nxl, nxr
417       IF ( psolver(1:7) == 'poisfft' )  THEN
418          DO  j = nys-1, nyn+1
419             DO  k = nzb, nzt+1
420                p(k,j,i) = tend(k,j,i)
421             ENDDO
422          ENDDO
423       ENDIF
424       DO  j = nys, nyn
425          DO  k = nzb_w_inner(j,i)+1, nzt
426             w(k,j,i) = w(k,j,i) - dt_3d * &
427                                   ( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1)
428          ENDDO
429          DO  k = nzb_u_inner(j,i)+1, nzt
430             u(k,j,i) = u(k,j,i) - dt_3d * ( tend(k,j,i) - tend(k,j,i-1) ) * ddx
431          ENDDO
432          DO  k = nzb_v_inner(j,i)+1, nzt
433             v(k,j,i) = v(k,j,i) - dt_3d * ( tend(k,j,i) - tend(k,j-1,i) ) * ddy
434          ENDDO
435
436!
437!--       Sum up the volume flow through the right and north boundary
438          IF ( conserve_volume_flow  .AND.  bc_lr == 'cyclic'  .AND. &
439               i == nx )  THEN
440             !$OMP CRITICAL
441             DO  k = nzb_2d(j,i) + 1, nzt
442                volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzu(k)
443             ENDDO
444             !$OMP END CRITICAL
445          ENDIF
446          IF ( conserve_volume_flow  .AND.  bc_ns == 'cyclic'  .AND. &
447               j == ny )  THEN
448             !$OMP CRITICAL
449             DO  k = nzb_2d(j,i) + 1, nzt
450                volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzu(k)
451             ENDDO
452             !$OMP END CRITICAL
453          ENDIF
454
455       ENDDO
456    ENDDO
457    !$OMP END PARALLEL
458
459!
460!-- Conserve the volume flow
461    IF ( conserve_volume_flow  .AND. &
462         ( bc_lr == 'cyclic'  .OR.  bc_ns == 'cyclic' ) )  THEN
463
464#if defined( __parallel )   
465       CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 2, MPI_REAL, &
466                           MPI_SUM, comm2d, ierr ) 
467#else
468       volume_flow = volume_flow_l 
469#endif   
470
471       volume_flow_offset = ( volume_flow_initial - volume_flow ) / &
472                            volume_flow_area
473
474       !$OMP PARALLEL PRIVATE (i,j,k)
475       !$OMP DO
476       DO  i = nxl, nxr
477          DO  j = nys, nyn
478             IF ( bc_lr == 'cyclic' )  THEN
479                DO  k = nzb_u_inner(j,i) + 1, nzt
480                   u(k,j,i) = u(k,j,i) + volume_flow_offset(1)
481                ENDDO
482             ENDIF
483             IF ( bc_ns == 'cyclic' )  THEN
484                DO  k = nzb_v_inner(j,i) + 1, nzt
485                   v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
486                ENDDO
487             ENDIF
488          ENDDO
489       ENDDO
490       !$OMP END PARALLEL
491
492    ENDIF
493
494!
495!-- Exchange of boundaries for the velocities
496    CALL exchange_horiz( u )
497    CALL exchange_horiz( v )
498    CALL exchange_horiz( w )
499
500!
501!-- Compute the divergence of the corrected velocity field,
502!-- a possible PE-sum is computed in flow_statistics
503    CALL cpu_log( log_point_s(1), 'divergence', 'start' )
504    sums_divnew_l = 0.0
505
506!
507!-- d must be reset to zero because it can contain nonzero values below the
508!-- topography
509    IF ( topography /= 'flat' )  d = 0.0
510
511    localsum  = 0.0
512    threadsum = 0.0
513
514    !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
515    !$OMP DO SCHEDULE( STATIC )
516#if defined( __ibm )
517    DO  i = nxl, nxr
518       DO  j = nys, nyn
519          DO  k = nzb_s_inner(j,i)+1, nzt
520             d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
521                        ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
522                        ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
523          ENDDO
524          DO  k = nzb+1, nzt
525             threadsum = threadsum + ABS( d(k,j,i) )
526          ENDDO
527       ENDDO
528    ENDDO
529#else
530    DO  i = nxl, nxr
531       DO  j = nys, nyn
532          DO  k = nzb_s_inner(j,i)+1, nzt
533             d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
534                        ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
535                        ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
536             threadsum = threadsum + ABS( d(k,j,i) )
537          ENDDO
538       ENDDO
539    ENDDO
540#endif
541    localsum = localsum + threadsum
542    !$OMP END PARALLEL
543
544!
545!-- For completeness, set the divergence sum of all statistic regions to those
546!-- of the total domain
547    sums_divnew_l(0:statistic_regions) = localsum
548
549    CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
550
551    CALL cpu_log( log_point(8), 'pres', 'stop' )
552
553
554 END SUBROUTINE pres
Note: See TracBrowser for help on using the repository browser.