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

Last change on this file since 1 was 1, checked in by raasch, 15 years ago

Initial repository layout and content

File size: 19.3 KB
Line 
1 SUBROUTINE pres
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Log: pres.f90,v $
11! Revision 1.25  2006/04/26 13:26:12  raasch
12! OpenMP optimization (+localsum, threadsum)
13!
14! Revision 1.24  2006/02/23 12:50:08  raasch
15! nzb replaced by nzb_.._inner, optional volume flow control, calculation of
16! divergence sum reordered due to optimization
17!
18! Revision 1.23  2005/05/18 15:53:44  raasch
19! call conditions for poisfft_hybrid modified
20!
21! Revision 1.22  2005/03/26 20:57:24  raasch
22! Adjustments for non-cyclic boundary conditions: arguments added to routine
23! exchange_horiz, mean horizontal wind parallel to the outflow is calculated.
24!
25! Revision 1.21  2004/04/30 12:45:00  raasch
26! Most of the arguments removed from poisfft call, +module poisfft_mod,
27! poisfft is also used for 1d-decompositions, more memory has to be
28! allocated for tend, when an enlarged domain is used for transposition
29!
30! Revision 1.20  2003/03/16 09:46:02  raasch
31! Two underscores (_) are placed in front of all define-strings
32!
33! Revision 1.19  2003/03/14 13:46:32  raasch
34! Optimization of loops (IBM-version and vectorizable version)
35!
36! Revision 1.18  2003/03/12 16:38:20  raasch
37! Simulated_time has been replaced by sums(nzb+1,4) in first if clause
38! in order to avoid run time errors due to compiler bugs on ibm and nec
39! machines, due to optimization, dividing d by dt_3d is now done in the
40! first loop (not in a seperated loop),
41! error removed: tend=p added after calling sor method
42!
43! Revision 1.17  2002/06/11 13:17:03  raasch
44! Hybrid solver included, array operations on d removed due to bad performance
45! on IBM, loops including the array d are combined, OpenMP directives added.
46! Total perturbation pressure is no more obtained iteratively (p=p+tend) from
47! the current perturbation pressure.
48!
49! Revision 1.16  2001/09/04  12:05:52  12:05:52  raasch (Siegfried Raasch)
50! Value of perturbation pressure is no more obtained iteratively (p = p + tend,
51! FFT-solver only)
52!
53! Revision 1.15  2001/07/20 13:12:54  raasch
54! Multigrid method included
55!
56! Revision 1.14  2001/03/30 07:47:29  raasch
57! All arguments removed, dp replaced by tend, array work removed from
58! sor argument list,
59! Translation of remaining German identifiers (variables, subroutines, etc.)
60!
61! Revision 1.13  2001/01/22 07:53:53  raasch
62! Module test_variables removed
63!
64! Revision 1.12  2000/01/26 13:23:32  letzel
65! All comments translated into English
66!
67! Revision 1.11  1998/09/22 17:28:14  raasch
68! dp war bei SOR-Aufrufen bisher nicht initialisiert, entsprechend geaendert
69!
70! Revision 1.10  1998/07/06 12:30:22  raasch
71! + USE test_variables
72!
73! Revision 1.9  1998/04/15 11:23:20  raasch
74! Testweise Zusatzbedingung fuer den Druck am unteren Rand bei inhomogener
75! Oberflaechentemperatur (pt wird uebergeben)
76!
77! Revision 1.8  1998/03/30 11:37:59  raasch
78! Divergenzsummenberechnung (PEs) nach flow_statistics verlagert
79!
80! Revision 1.7  1998/03/25 13:56:16  raasch
81! dt in dt_3d umbenannt
82!
83! Revision 1.6  1997/09/16 06:38:37  raasch
84! Ueberfluessige Initialisierung von d entfernt
85!
86! Revision 1.5  1997/09/12 06:29:56  raasch
87! Randbedingungen umbenannt
88!
89! Revision 1.4  1997/09/09 08:30:27  raasch
90! Kehrwerte der Gitterweiten implementiert
91!
92! Revision 1.3  1997/08/26 06:34:16  raasch
93! Gesamtstoerdruck p wird in Bew.gleichungen mitgerechnet und ergibt sich
94! durch Summation des aktuellen Stoerdrucks dp, mit dem hier korrigiert wird
95!
96! Revision 1.2  1997/08/11 06:25:24  raasch
97! Felder werden ueber Formalparameterliste uebergeben.
98!
99! Revision 1.1  1997/07/24 11:24:44  raasch
100! Initial revision
101!
102!
103! Description:
104! ------------
105! Compute the divergence of the provisional velocity field. Solve the Poisson
106! equation for the perturbation pressure. Compute the final velocities using
107! this perturbation pressure. Compute the remaining divergence.
108!------------------------------------------------------------------------------!
109
110    USE arrays_3d
111    USE constants
112    USE control_parameters
113    USE cpulog
114    USE grid_variables
115    USE indices
116    USE interfaces
117    USE pegrid
118    USE poisfft_mod
119    USE poisfft_hybrid_mod
120    USE statistics
121
122    IMPLICIT NONE
123
124    INTEGER ::  i, j, k, sr
125
126    REAL    ::  localsum, threadsum
127
128    REAL, DIMENSION(1:2) ::  volume_flow_l, volume_flow_offset
129
130
131    CALL cpu_log( log_point(8), 'pres', 'start' )
132
133!
134!-- Multigrid method needs additional grid points for the divergence array
135    IF ( psolver == 'multigrid' )  THEN
136       DEALLOCATE( d )
137       ALLOCATE( d(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
138    ENDIF
139
140!
141!-- Compute the divergence of the provisional velocity field.
142    CALL cpu_log( log_point_s(1), 'divergence', 'start' )
143
144    IF ( psolver == 'multigrid' )  THEN
145       !$OMP PARALLEL DO SCHEDULE( STATIC )
146       DO  i = nxl-1, nxr+1
147          DO  j = nys-1, nyn+1
148             DO  k = nzb, nzt+1
149                d(k,j,i) = 0.0
150             ENDDO
151          ENDDO
152       ENDDO
153    ELSE
154       !$OMP PARALLEL DO SCHEDULE( STATIC )
155       DO  i = nxl, nxra
156          DO  j = nys, nyna
157             DO  k = nzb+1, nzta
158                d(k,j,i) = 0.0
159             ENDDO
160          ENDDO
161       ENDDO
162    ENDIF
163
164    localsum  = 0.0
165    threadsum = 0.0
166
167#if defined( __ibm )
168    !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
169    !$OMP DO SCHEDULE( STATIC )
170    DO  i = nxl, nxr
171       DO  j = nys, nyn
172          DO  k = nzb_s_inner(j,i)+1, nzt
173             d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
174                        ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
175                        ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
176          ENDDO
177!
178!--       Additional pressure boundary condition at the bottom boundary for
179!--       inhomogeneous Prandtl layer heat fluxes and temperatures, respectively
180!--       dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0.
181!--       This condition must not be applied at the start of a run, because then
182!--       flow_statistics has not yet been called and thus sums = 0.
183          IF ( ibc_p_b == 2  .AND.  sums(nzb+1,4) /= 0.0 )  THEN
184             k = nzb_s_inner(j,i)
185             d(k+1,j,i) = d(k+1,j,i) + (                                     &
186                                         ( usws(j,i+1) - usws(j,i) ) * ddx   &
187                                       + ( vsws(j+1,i) - vsws(j,i) ) * ddy   &
188                                       - g * ( pt(k+1,j,i) - sums(k+1,4) ) / &
189                                         sums(k+1,4)                         &
190                                       ) * ddzw(k+1)
191          ENDIF
192
193!
194!--       Compute possible PE-sum of divergences for flow_statistics
195          DO  k = nzb_s_inner(j,i)+1, nzt
196             threadsum = threadsum + ABS( d(k,j,i) )
197          ENDDO
198
199!
200!--       Velocity corrections are made with Euler step size. Right hand side
201!--       of Poisson equation has to be set appropriately
202          DO  k = nzb_s_inner(j,i)+1, nzt
203             d(k,j,i) = d(k,j,i) / dt_3d
204          ENDDO
205
206       ENDDO
207    ENDDO
208
209    localsum = localsum + threadsum
210    !$OMP END PARALLEL
211#else
212    IF ( ibc_p_b == 2 .AND. sums(nzb+1,4) /= 0.0 )  THEN
213       !$OMP PARALLEL PRIVATE (i,j,k)
214       !$OMP DO SCHEDULE( STATIC )
215       DO  i = nxl, nxr
216          DO  j = nys, nyn
217             DO  k = nzb_s_inner(j,i)+1, nzt
218                d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
219                           ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
220                           ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
221             ENDDO
222          ENDDO
223!
224!--       Additional pressure boundary condition at the bottom boundary for
225!--       inhomogeneous Prandtl layer heat fluxes and temperatures, respectively
226!--       dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0.
227!--       This condition must not be applied at the start of a run, because then
228!--       flow_statistics has not yet been called and thus sums = 0.
229          DO  j = nys, nyn
230              k = nzb_s_inner(j,i)
231              d(k+1,j,i) = d(k+1,j,i) + (                        &
232                             ( usws(j,i+1) - usws(j,i) ) * ddx   &
233                           + ( vsws(j+1,i) - vsws(j,i) ) * ddy   &
234                           - g * ( pt(k+1,j,i) - sums(k+1,4) ) / &
235                             sums(k+1,4)                         &
236                                        ) * ddzw(k+1)
237          ENDDO
238       ENDDO
239       !$OMP END PARALLEL
240
241    ELSE
242
243       !$OMP PARALLEL PRIVATE (i,j,k)
244       !$OMP DO SCHEDULE( STATIC )
245       DO  i = nxl, nxr
246          DO  j = nys, nyn
247             DO  k = nzb_s_inner(j,i)+1, nzt
248                d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
249                           ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
250                           ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
251             ENDDO
252          ENDDO
253       ENDDO
254       !$OMP END PARALLEL
255
256    ENDIF
257
258!
259!-- Compute possible PE-sum of divergences for flow_statistics
260    !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
261    !$OMP DO SCHEDULE( STATIC )
262    DO  i = nxl, nxr
263       DO  j = nys, nyn
264          DO  k = nzb+1, nzt
265             threadsum = threadsum + ABS( d(k,j,i) )
266          ENDDO
267       ENDDO
268    ENDDO
269    localsum = localsum + threadsum
270    !$OMP END PARALLEL
271
272!
273!-- Velocity corrections are made with Euler step size. Right hand side
274!-- of Poisson equation has to be set appropriately
275    !$OMP DO SCHEDULE( STATIC )
276    DO  i = nxl, nxr
277       DO  j = nys, nyn
278          DO  k = nzb_s_inner(j,i)+1, nzt
279             d(k,j,i) = d(k,j,i) / dt_3d
280          ENDDO
281       ENDDO
282    ENDDO
283#endif
284
285!
286!-- For completeness, set the divergence sum of all statistic regions to those
287!-- of the total domain
288    sums_divold_l(0:statistic_regions) = localsum
289
290!
291!-- Determine absolute minimum/maximum (only for test cases, therefore as
292!-- comment line)
293!    CALL global_min_max( nzb+1, nzt, nys, nyn, nxl, nxr, d, 'abs', divmax, &
294!                         divmax_ijk )
295
296    CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
297
298!
299!-- Compute the pressure perturbation solving the Poisson equation
300    IF ( psolver(1:7) == 'poisfft' )  THEN
301
302!
303!--    Enlarge the size of tend, used as a working array for the transpositions
304       IF ( nxra > nxr  .OR.  nyna > nyn  .OR.  nza > nz )  THEN
305          DEALLOCATE( tend )
306          ALLOCATE( tend(1:nza,nys:nyna,nxl:nxra) )
307       ENDIF
308
309!
310!--    Solve Poisson equation via FFT and solution of tridiagonal matrices
311       IF ( psolver == 'poisfft' )  THEN
312!
313!--       Solver for 2d-decomposition
314          CALL poisfft( d, tend )
315       ELSEIF ( psolver == 'poisfft_hybrid' )  THEN 
316!
317!--       Solver for 1d-decomposition (using MPI and OpenMP).
318!--       The old hybrid-solver is still included here, as long as there
319!--       are some optimization problems in poisfft
320          CALL poisfft_hybrid( d )
321       ENDIF
322
323!
324!--    Resize tend to its normal size
325       IF ( nxra > nxr  .OR.  nyna > nyn  .OR.  nza > nz )  THEN
326          DEALLOCATE( tend )
327          ALLOCATE( tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
328       ENDIF
329
330!
331!--    Store computed perturbation pressure and set boundary condition in
332!--    z-direction
333       !$OMP PARALLEL DO
334       DO  i = nxl, nxr
335          DO  j = nys, nyn
336             DO  k = nzb+1, nzt
337                tend(k,j,i) = d(k,j,i)
338             ENDDO
339          ENDDO
340       ENDDO
341
342!
343!--    Bottom boundary:
344!--    This condition is only required for internal output. The pressure
345!--    gradient (dp(nzb+1)-dp(nzb))/dz is not used anywhere else.
346       IF ( ibc_p_b == 1 )  THEN
347!
348!--       Neumann (dp/dz = 0)
349          !$OMP PARALLEL DO
350          DO  i = nxl-1, nxr+1
351             DO  j = nys-1, nyn+1
352                tend(nzb_s_inner(j,i),j,i) = tend(nzb_s_inner(j,i)+1,j,i)
353             ENDDO
354          ENDDO
355
356       ELSEIF ( ibc_p_b == 2 )  THEN
357!
358!--       Neumann condition for inhomogeneous surfaces,
359!--       here currently still in the form of a zero gradient. Actually
360!--       dp/dz = -(dtau13/dx + dtau23/dy) + g*pt'/pt0 would have to be used for
361!--       the computation (cf. above: computation of divergences).
362          !$OMP PARALLEL DO
363          DO  i = nxl-1, nxr+1
364             DO  j = nys-1, nyn+1
365                tend(nzb_s_inner(j,i),j,i) = tend(nzb_s_inner(j,i)+1,j,i)
366             ENDDO
367          ENDDO
368
369       ELSE
370!
371!--       Dirichlet
372          !$OMP PARALLEL DO
373          DO  i = nxl-1, nxr+1
374             DO  j = nys-1, nyn+1
375                tend(nzb_s_inner(j,i),j,i) = 0.0
376             ENDDO
377          ENDDO
378
379       ENDIF
380
381!
382!--    Top boundary
383       IF ( ibc_p_t == 1 )  THEN
384!
385!--       Neumann
386          !$OMP PARALLEL DO
387          DO  i = nxl-1, nxr+1
388             DO  j = nys-1, nyn+1
389                tend(nzt+1,j,i) = tend(nzt,j,i)
390             ENDDO
391          ENDDO
392
393       ELSE
394!
395!--       Dirichlet
396          !$OMP PARALLEL DO
397          DO  i = nxl-1, nxr+1
398             DO  j = nys-1, nyn+1
399                tend(nzt+1,j,i) = 0.0
400             ENDDO
401          ENDDO
402
403       ENDIF
404
405!
406!--    Exchange boundaries for p
407       CALL exchange_horiz( tend, 0, 0 )
408     
409    ELSEIF ( psolver == 'sor' )  THEN
410
411!
412!--    Solve Poisson equation for perturbation pressure using SOR-Red/Black
413!--    scheme
414       CALL sor( d, ddzu, ddzw, p )
415       tend = p
416
417    ELSEIF ( psolver == 'multigrid' )  THEN
418
419!
420!--    Solve Poisson equation for perturbation pressure using Multigrid scheme,
421!--    array tend is used to store the residuals
422       CALL poismg( tend )
423 
424!
425!--    Restore perturbation pressure on tend because this array is used
426!--    further below to correct the velocity fields
427       tend = p
428
429    ENDIF
430
431!
432!-- Store perturbation pressure on array p, used in the momentum equations
433    IF ( psolver(1:7) == 'poisfft' )  THEN
434!
435!--    Here, only the values from the left and right boundaries are copied
436!--    The remaining values are copied in the following loop due to speed
437!--    optimization
438       !$OMP PARALLEL DO
439       DO  j = nys-1, nyn+1
440          DO  k = nzb, nzt+1
441             p(k,j,nxl-1) = tend(k,j,nxl-1)
442             p(k,j,nxr+1) = tend(k,j,nxr+1)
443          ENDDO
444       ENDDO
445    ENDIF
446
447!
448!-- Correction of the provisional velocities with the current perturbation
449!-- pressure just computed
450    IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  uvmean_outflow_l = 0.0
451    IF ( conserve_volume_flow )  THEN
452       volume_flow_l(1) = 0.0
453       volume_flow_l(2) = 0.0
454    ENDIF
455    !$OMP PARALLEL PRIVATE (i,j,k)
456    !$OMP DO
457    DO  i = nxl, nxr
458       IF ( psolver(1:7) == 'poisfft' )  THEN
459          DO  j = nys-1, nyn+1
460             DO  k = nzb, nzt+1
461                p(k,j,i) = tend(k,j,i)
462             ENDDO
463          ENDDO
464       ENDIF
465       DO  j = nys, nyn
466          DO  k = nzb_w_inner(j,i)+1, nzt
467             w(k,j,i) = w(k,j,i) - dt_3d * &
468                                   ( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1)
469          ENDDO
470          DO  k = nzb_u_inner(j,i)+1, nzt
471             u(k,j,i) = u(k,j,i) - dt_3d * ( tend(k,j,i) - tend(k,j,i-1) ) * ddx
472          ENDDO
473          DO  k = nzb_v_inner(j,i)+1, nzt
474             v(k,j,i) = v(k,j,i) - dt_3d * ( tend(k,j,i) - tend(k,j-1,i) ) * ddy
475          ENDDO
476
477!
478!--       Sum up the horizontal velocity along the outflow plane (in case
479!--       of non-cyclic boundary conditions). The respective mean velocity
480!--       is calculated from this in routine boundary_conds.
481          IF ( outflow_l  .AND.  i == nxl )  THEN
482             !$OMP CRITICAL
483             DO  k = nzb, nzt+1
484                uvmean_outflow_l(k) = uvmean_outflow_l(k) + v(k,j,nxl)
485             ENDDO
486             !$OMP END CRITICAL
487          ELSEIF ( outflow_r  .AND.  i == nxr )  THEN
488             !$OMP CRITICAL
489             DO  k = nzb, nzt+1
490                uvmean_outflow_l(k) = uvmean_outflow_l(k) + v(k,j,nxr)
491             ENDDO
492             !$OMP END CRITICAL
493          ELSEIF ( outflow_s  .AND.  j == nys )  THEN
494             !$OMP CRITICAL
495             DO  k = nzb, nzt+1
496                uvmean_outflow_l(k) = uvmean_outflow_l(k) + u(k,nys,i)
497             ENDDO
498             !$OMP END CRITICAL
499          ELSEIF ( outflow_n  .AND.  j == nyn )  THEN
500             !$OMP CRITICAL
501             DO  k = nzb, nzt+1
502                uvmean_outflow_l(k) = uvmean_outflow_l(k) + u(k,nyn,i)
503             ENDDO
504             !$OMP END CRITICAL
505          ENDIF
506
507!
508!--       Sum up the volume flow through the right and north boundary
509          IF ( conserve_volume_flow  .AND.  i == nx )  THEN
510             !$OMP CRITICAL
511             DO  k = nzb_2d(j,i) + 1, nzt
512                volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzu(k)
513             ENDDO
514             !$OMP END CRITICAL
515          ENDIF
516          IF ( conserve_volume_flow  .AND.  j == ny )  THEN
517             !$OMP CRITICAL
518             DO  k = nzb_2d(j,i) + 1, nzt
519                volume_flow_l(2) = volume_flow_l(2) + v(k,j,i) * dzu(k)
520             ENDDO
521             !$OMP END CRITICAL
522          ENDIF
523
524       ENDDO
525    ENDDO
526    !$OMP END PARALLEL
527
528!
529!-- Conserve the volume flow
530    IF ( conserve_volume_flow )  THEN
531
532#if defined( __parallel )   
533       CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 2, MPI_REAL, &
534                           MPI_SUM, comm2d, ierr ) 
535#else
536       volume_flow = volume_flow_l 
537#endif   
538
539       volume_flow_offset = ( volume_flow_initial - volume_flow ) / &
540                            volume_flow_area
541
542       !$OMP PARALLEL PRIVATE (i,j,k)
543       !$OMP DO
544       DO  i = nxl, nxr
545          DO  j = nys, nyn
546             DO  k = nzb_u_inner(j,i) + 1, nzt
547                u(k,j,i) = u(k,j,i) + volume_flow_offset(1)
548             ENDDO
549             DO  k = nzb_v_inner(j,i) + 1, nzt
550                v(k,j,i) = v(k,j,i) + volume_flow_offset(2)
551             ENDDO
552          ENDDO
553       ENDDO
554       !$OMP END PARALLEL
555
556    ENDIF
557
558!
559!-- Exchange of boundaries for the velocities
560    CALL exchange_horiz( u, uxrp,    0 )
561    CALL exchange_horiz( v,    0, vynp )
562    CALL exchange_horiz( w,    0,    0 )
563
564!
565!-- Compute the divergence of the corrected velocity field,
566!-- a possible PE-sum is computed in flow_statistics
567    CALL cpu_log( log_point_s(1), 'divergence', 'start' )
568    sums_divnew_l = 0.0
569
570!
571!-- d must be reset to zero because it can contain nonzero values below the
572!-- topography
573    IF ( topography /= 'flat' )  d = 0.0
574
575    localsum  = 0.0
576    threadsum = 0.0
577
578    !$OMP PARALLEL PRIVATE (i,j,k) FIRSTPRIVATE(threadsum) REDUCTION(+:localsum)
579    !$OMP DO SCHEDULE( STATIC )
580#if defined( __ibm )
581    DO  i = nxl, nxr
582       DO  j = nys, nyn
583          DO  k = nzb_s_inner(j,i)+1, nzt
584             d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
585                        ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
586                        ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
587          ENDDO
588          DO  k = nzb+1, nzt
589             threadsum = threadsum + ABS( d(k,j,i) )
590          ENDDO
591       ENDDO
592    ENDDO
593#else
594    DO  i = nxl, nxr
595       DO  j = nys, nyn
596          DO  k = nzb_s_inner(j,i)+1, nzt
597             d(k,j,i) = ( u(k,j,i+1) - u(k,j,i) ) * ddx + &
598                        ( v(k,j+1,i) - v(k,j,i) ) * ddy + &
599                        ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
600             threadsum = threadsum + ABS( d(k,j,i) )
601          ENDDO
602       ENDDO
603    ENDDO
604#endif
605    localsum = localsum + threadsum
606    !$OMP END PARALLEL
607
608!
609!-- For completeness, set the divergence sum of all statistic regions to those
610!-- of the total domain
611    sums_divnew_l(0:statistic_regions) = localsum
612
613    CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
614
615    CALL cpu_log( log_point(8), 'pres', 'stop' )
616
617
618 END SUBROUTINE pres
Note: See TracBrowser for help on using the repository browser.