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

Last change on this file since 622 was 622, checked in by raasch, 11 years ago

New:
---

Optional barriers included in order to speed up collective operations
MPI_ALLTOALL and MPI_ALLREDUCE. This feature is controlled with new initial
parameter collective_wait. Default is .FALSE, but .TRUE. on SGI-type
systems. (advec_particles, advec_s_bc, buoyancy, check_for_restart,
cpu_statistics, data_output_2d, data_output_ptseries, flow_statistics,
global_min_max, inflow_turbulence, init_3d_model, init_particles, init_pegrid,
init_slope, parin, pres, poismg, set_particle_attributes, timestep,
read_var_list, user_statistics, write_compressed, write_var_list)

Adjustments for Kyushu Univ. (lcrte, ibmku). Concerning hybrid
(MPI/openMP) runs, the number of openMP threads per MPI tasks can now
be given as an argument to mrun-option -O. (mbuild, mrun, subjob)

Changed:


Initialization of the module command changed for SGI-ICE/lcsgi (mbuild, subjob)

Errors:


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