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

Last change on this file since 76 was 76, checked in by raasch, 14 years ago

changes for Neumann p conditions both at bottom and top

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