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

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

New:
---
Allows runs for a coupled atmosphere-ocean LES,
coupling frequency is controlled by new d3par-parameter dt_coupling,
the coupling mode (atmosphere_to_ocean or ocean_to_atmosphere) for the
respective processes is read from environment variable coupling_mode,
which is set by the mpiexec-command,
communication between the two models is done using the intercommunicator
comm_inter,
local files opened by the ocean model get the additional suffic "_O".
Assume saturation at k=nzb_s_inner(j,i) for atmosphere coupled to ocean.

A momentum flux can be set as top boundary condition using the new
inipar parameter top_momentumflux_u|v.

Non-cyclic boundary conditions can be used along all horizontal directions.

Quantities w*p* and w"e can be output as vertical profiles.

Initial profiles are reset to constant profiles in case that initializing_actions /= 'set_constant_profiles'. (init_rankine)

Optionally calculate km and kh from initial TKE e_init.

Changed:


Remaining variables iran changed to iran_part (advec_particles, init_particles).

In case that the presure solver is not called for every Runge-Kutta substep
(call_psolver_at_all_substeps = .F.), it is called after the first substep
instead of the last. In that case, random perturbations are also added to the
velocity field after the first substep.

Initialization of km,kh = 0.00001 for ocean = .T. (for ocean = .F. it remains 0.01).

Allow data_output_pr= q, wq, w"q", w*q* for humidity = .T. (instead of cloud_physics = .T.).

Errors:


Bugs from code parts for non-cyclic boundary conditions are removed: loops for
u and v are starting from index nxlu, nysv, respectively. The radiation boundary
condition is used for every Runge-Kutta substep. Velocity phase speeds for
the radiation boundary conditions are calculated for the first Runge-Kutta
substep only and reused for the further substeps. New arrays c_u, c_v, and c_w
are defined for this purpose. Several index errors are removed from the
radiation boundary condition code parts. Upper bounds for calculating
u_0 and v_0 (in production_e) are nxr+1 and nyn+1 because otherwise these
values are not available in case of non-cyclic boundary conditions.

+dots_num_palm in module user, +module netcdf_control in user_init (both in user_interface)

Bugfix: wrong sign removed from the buoyancy production term in the case use_reference = .T. (production_e)

Bugfix: Error message concerning output of particle concentration (pc) modified (check_parameters).

Bugfix: Rayleigh damping for ocean fixed.

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