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

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

preliminary update of bugfixes and extensions for non-cyclic BCs

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