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

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

last commit documented

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