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

Last change on this file since 169 was 151, checked in by raasch, 17 years ago

preliminary update for the turbulence recycling method

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