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

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

openmp bugfixes found by NEC benchmarker

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