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

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

typo in file headers removed

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