source: palm/tags/release-3.2/SOURCE/pres.f90 @ 2238

Last change on this file since 2238 was 77, checked in by raasch, 17 years ago

New:
---

particle reflection from vertical walls implemented, particle SGS model adjusted to walls

Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.

new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files

new d3par-parameter dt_max to define the maximum value for the allowed timestep

new inipar-parameter loop_optimization to control the loop optimization method

new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).

new user interface user_advec_particles

new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays

topography height informations are stored on arrays zu_s_inner and zw_w_inner and output to the 2d/3d NetCDF files

samples added to the user interface which show how to add user-define time series quantities.

calculation/output of precipitation amount, precipitation rate and z0 (by setting "pra*", "prr*", "z0*" with data_output). The time interval on which the precipitation amount is defined is set by new d3par-parameter precipitation_amount_interval

unit 9 opened for debug output (file DEBUG_<pe#>)

Makefile, advec_particles, average_3d_data, buoyancy, calc_precipitation, check_open, check_parameters, data_output_2d, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, impact_of_latent_heat, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, read_3d_binary, sum_up_3d_data, user_interface, write_var_list, write_3d_binary

New: wall_fluxes

Changed:


General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, volume flow control is added for the outflow boundary (currently only for the north boundary!!), additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d

The FFT-method for solving the Poisson-equation is now working with Neumann boundary conditions both at the bottom and the top. This requires adjustments of the tridiagonal coefficients and subtracting the horizontally averaged mean from the vertical velocity field.

+age_m in particle_type

Particles-package is now part of the default code ("-p particles" is not needed any more).

Move call of user_actions( 'after_integration' ) below increment of times
and counters. user_actions is now called for each statistic region and has as an argument the number of the respective region (sr)

d3par-parameter data_output_ts removed. Timeseries output for "profil" removed. Timeseries are now switched on by dt_dots. Timeseries data is collected in flow_statistics.

Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.

q is not allowed to become negative (prognostic_equations).

poisfft_init is only called if fft-solver is switched on (init_pegrid).

d3par-parameter moisture renamed to humidity.

Subversion global revision number is read from mrun and added to the run description header and to the run control (_rc) file.

vtk directives removed from main program.

The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.

advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poisfft, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary

Errors:


Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model

Bugfix in sample for reading user defined data from restart file (user_init)

Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)

Check for possible negative humidities in the initial humidity profile.

in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
# case of .mod files

Makefile
check_parameters, init_1d_model, init_3d_model, user_interface

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