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

Last change on this file since 44 was 4, checked in by raasch, 18 years ago

Id keyword set as property for all *.f90 files

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