Index: palm/trunk/SOURCE/pres.f90
===================================================================
--- palm/trunk/SOURCE/pres.f90	(revision 695)
+++ palm/trunk/SOURCE/pres.f90	(revision 707)
@@ -4,5 +4,10 @@
 ! Current revisions:
 ! -----------------
-! 
+! Calculation of weighted average of p is now handled in the same way
+! regardless of the number of ghost layers (advection scheme),
+! multigrid and sor method are using p_loc for iterative advancements of
+! pressure, 
+! localsum calculation modified for proper OpenMP reduction,
+! bc_lr/ns replaced by bc_lr/ns_cyc
 !
 ! Former revisions:
@@ -99,20 +104,33 @@
 
 !
-!-- Multigrid method expects 1 additional grid point for the arrays 
-!-- d, tend and p
+!-- Multigrid method expects array d to have one ghost layer.
+!-- 
     IF ( psolver == 'multigrid' )  THEN
      
        DEALLOCATE( d )
        ALLOCATE( d(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ) 
+
+!
+!--    Since p is later used to hold the weighted average of the substeps, it
+!--    cannot be used in the iterative solver. Therefore, its initial value is
+!--    stored on p_loc, which is then iteratively advanced in every substep.
+       IF ( intermediate_timestep_count == 1 )  THEN
+          DO  i = nxl-1, nxr+1
+             DO  j = nys-1, nyn+1
+                DO  k = nzb, nzt+1
+                   p_loc(k,j,i) = p(k,j,i)
+                ENDDO
+             ENDDO
+          ENDDO
+       ENDIF
        
-       IF ( ws_scheme_mom  .OR. ws_scheme_sca )  THEN
-       
-          DEALLOCATE( tend )
-          ALLOCATE( tend(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
-          DEALLOCATE( p )
-          ALLOCATE( p(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
-          
-       ENDIF
-       
+    ELSEIF ( psolver == 'sor'  .AND.  intermediate_timestep_count == 1 )  THEN
+
+!
+!--    Since p is later used to hold the weighted average of the substeps, it
+!--    cannot be used in the iterative solver. Therefore, its initial value is
+!--    stored on p_loc, which is then iteratively advanced in every substep.
+       p_loc = p
+
     ENDIF
 
@@ -297,6 +315,6 @@
     ENDDO
 
-    localsum = ( localsum + threadsum ) * dt_3d * &
-               weight_pres(intermediate_timestep_count)
+    localsum = localsum + threadsum * dt_3d * &
+                          weight_pres(intermediate_timestep_count)
 
     !$OMP END PARALLEL
@@ -362,6 +380,6 @@
        ENDDO
     ENDDO
-    localsum = ( localsum + threadsum ) * dt_3d                               &
-               * weight_pres(intermediate_timestep_count)
+    localsum = localsum + threadsum * dt_3d * &
+                          weight_pres(intermediate_timestep_count)
     !$OMP END PARALLEL
 #endif
@@ -371,10 +389,4 @@
 !-- of the total domain
     sums_divold_l(0:statistic_regions) = localsum
-
-!
-!-- Determine absolute minimum/maximum (only for test cases, therefore as
-!-- comment line)
-!    CALL global_min_max( nzb+1, nzt, nys, nyn, nxl, nxr, d, 'abs', divmax, &
-!                        divmax_ijk )
 
     CALL cpu_log( log_point_s(1), 'divergence', 'stop' )
@@ -496,6 +508,6 @@
 !--    Solve Poisson equation for perturbation pressure using SOR-Red/Black
 !--    scheme
-       CALL sor( d, ddzu_pres, ddzw, p )
-       tend = p
+       CALL sor( d, ddzu_pres, ddzw, p_loc )
+       tend = p_loc
 
     ELSEIF ( psolver == 'multigrid' )  THEN
@@ -506,64 +518,62 @@
 !--    to discern data exchange in multigrid ( 1 ghostpoint ) and normal grid
 !--    ( nbgp ghost points ).
-       exchange_mg = .TRUE.
        CALL poismg( tend )
-       exchange_mg = .FALSE.
+
 !
 !--    Restore perturbation pressure on tend because this array is used
 !--    further below to correct the velocity fields
-
-       tend = p
-       IF( ws_scheme_mom .OR. ws_scheme_sca )  THEN
-!       
-!--       Allocate p to its normal size and restore pressure.         
-          DEALLOCATE( p )
-          ALLOCATE( p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-          
-       ENDIF
-
-    ENDIF
-
-!
-!-- Store perturbation pressure on array p, used in the momentum equations
-    IF ( psolver(1:7) == 'poisfft' )  THEN
-
-       IF ( intermediate_timestep_count == 1 )  THEN
-          !$OMP PARALLEL PRIVATE (i,j,k)
-          !$OMP DO
-          DO  i = nxlg, nxrg
-             DO  j = nysg, nyng
-                DO  k = nzb, nzt+1
-                   p(k,j,i) = tend(k,j,i)                                     &
-                            * weight_substep(intermediate_timestep_count)
-                ENDDO
-             ENDDO
-          ENDDO
-          !$OMP END PARALLEL
-          
-       ELSE 
-          !$OMP PARALLEL PRIVATE (i,j,k)
-          !$OMP DO
-          DO  i = nxlg, nxrg
-             DO  j = nysg, nyng
-                DO  k = nzb, nzt+1
-                   p(k,j,i) = p(k,j,i) + tend(k,j,i)                          &
-                            * weight_substep(intermediate_timestep_count)
-                ENDDO
-             ENDDO
-          ENDDO
-          !$OMP END PARALLEL
-          
-       ENDIF
+       DO  i = nxl-1, nxr+1
+          DO  j = nys-1, nyn+1
+             DO  k = nzb, nzt+1
+                tend(k,j,i) = p_loc(k,j,i)
+             ENDDO
+          ENDDO
+       ENDDO
+
+    ENDIF
+
+!
+!-- Store perturbation pressure on array p, used for pressure data output.
+!-- Ghost layers are added in the output routines (except sor-method: see below)
+    IF ( intermediate_timestep_count == 1 )  THEN
+       !$OMP PARALLEL PRIVATE (i,j,k)
+       !$OMP DO
+       DO  i = nxl-1, nxr+1
+          DO  j = nys-1, nyn+1
+             DO  k = nzb, nzt+1
+                p(k,j,i) = tend(k,j,i) * &
+                           weight_substep(intermediate_timestep_count)
+             ENDDO
+          ENDDO
+       ENDDO
+       !$OMP END PARALLEL
+
+    ELSE 
+       !$OMP PARALLEL PRIVATE (i,j,k)
+       !$OMP DO
+       DO  i = nxl-1, nxr+1
+          DO  j = nys-1, nyn+1
+             DO  k = nzb, nzt+1
+                p(k,j,i) = p(k,j,i) + tend(k,j,i) * &
+                           weight_substep(intermediate_timestep_count)
+             ENDDO
+          ENDDO
+       ENDDO
+       !$OMP END PARALLEL
+
+    ENDIF
        
-    ENDIF
+!
+!-- SOR-method needs ghost layers for the next timestep
+    IF ( psolver == 'sor' )  CALL exchange_horiz( p, nbgp )
 
 !
 !-- Correction of the provisional velocities with the current perturbation 
 !-- pressure just computed
-    IF ( conserve_volume_flow  .AND. &
-         ( bc_lr == 'cyclic'  .OR.  bc_ns == 'cyclic' ) )  THEN
+    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc  .OR.  bc_ns_cyc ) )  THEN
        volume_flow_l(1) = 0.0
        volume_flow_l(2) = 0.0
     ENDIF
+
     !$OMP PARALLEL PRIVATE (i,j,k)
     !$OMP DO
@@ -571,22 +581,22 @@
        DO  j = nys, nyn
           DO  k = nzb_w_inner(j,i)+1, nzt
-             w(k,j,i) = w(k,j,i) - dt_3d *                                    &
-                        ( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1)           &
-                        * weight_pres(intermediate_timestep_count)
+             w(k,j,i) = w(k,j,i) - dt_3d *                                 &
+                        ( tend(k+1,j,i) - tend(k,j,i) ) * ddzu(k+1) *      &
+                        weight_pres(intermediate_timestep_count)
           ENDDO
           DO  k = nzb_u_inner(j,i)+1, nzt
              u(k,j,i) = u(k,j,i) - dt_3d *                                 &
-                        ( tend(k,j,i) - tend(k,j,i-1) ) * ddx              &
-                        * weight_pres(intermediate_timestep_count)
+                        ( tend(k,j,i) - tend(k,j,i-1) ) * ddx *            &
+                        weight_pres(intermediate_timestep_count)
           ENDDO
           DO  k = nzb_v_inner(j,i)+1, nzt
              v(k,j,i) = v(k,j,i) - dt_3d *                                 &
-                        ( tend(k,j,i) - tend(k,j-1,i) ) * ddy              &
-                        * weight_pres(intermediate_timestep_count)
+                        ( tend(k,j,i) - tend(k,j-1,i) ) * ddy *            &
+                        weight_pres(intermediate_timestep_count)
           ENDDO                                                         
 !
 !--       Sum up the volume flow through the right and north boundary
-          IF ( conserve_volume_flow  .AND.  bc_lr == 'cyclic'  .AND. &
-               bc_ns == 'cyclic' .AND. i == nx )  THEN
+          IF ( conserve_volume_flow  .AND.  bc_lr_cyc  .AND.  bc_ns_cyc  .AND. &
+               i == nx )  THEN
              !$OMP CRITICAL
              DO  k = nzb_2d(j,i) + 1, nzt
@@ -595,6 +605,6 @@
              !$OMP END CRITICAL
           ENDIF
-          IF ( conserve_volume_flow  .AND.  bc_ns == 'cyclic'  .AND. &
-               bc_lr == 'cyclic' .AND. j == ny )  THEN
+          IF ( conserve_volume_flow  .AND.  bc_ns_cyc  .AND.  bc_lr_cyc  .AND. &
+               j == ny )  THEN
              !$OMP CRITICAL
              DO  k = nzb_2d(j,i) + 1, nzt
@@ -608,59 +618,7 @@
     !$OMP END PARALLEL
     
-    IF ( psolver == 'multigrid' .OR. psolver == 'sor' )  THEN       
-       IF ( intermediate_timestep_count == 1 .OR. simulated_time == 0)  THEN 
-          !$OMP PARALLEL PRIVATE (i,j,k)
-          !$OMP DO
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1
-                   p_sub(k,j,i) = tend(k,j,i)                                 &
-                                * weight_substep(intermediate_timestep_count)
-                ENDDO
-             ENDDO
-          ENDDO
-          !$OMP END PARALLEL
-       ELSE
-          !$OMP PARALLEL PRIVATE (i,j,k)
-          !$OMP DO
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1
-                   p_sub(k,j,i) = p_sub(k,j,i) + tend(k,j,i)                  &
-                                * weight_substep(intermediate_timestep_count)
-                ENDDO
-             ENDDO
-          ENDDO
-          !$OMP END PARALLEL
-       ENDIF
-          
-       IF ( intermediate_timestep_count == intermediate_timestep_count_max )  &
-          THEN
-          !$OMP PARALLEL PRIVATE (i,j,k)
-          !$OMP DO
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO  k = nzb, nzt+1
-                   p(k,j,i) = p_sub(k,j,i) 
-                ENDDO
-             ENDDO
-          ENDDO
-          !$OMP END PARALLEL
-       ENDIF
-    ENDIF
- 
-!
-!-- Resize tend to its normal size in case of multigrid and ws-scheme.
-    IF ( psolver == 'multigrid' .AND. ( ws_scheme_mom        &
-                                   .OR. ws_scheme_sca ) )  THEN
-       DEALLOCATE( tend )
-       ALLOCATE( tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-    ENDIF
-
-
 !
 !-- Conserve the volume flow
-    IF ( conserve_volume_flow  .AND. &
-         ( bc_lr == 'cyclic'  .AND.  bc_ns == 'cyclic' ) )  THEN
+    IF ( conserve_volume_flow  .AND.  ( bc_lr_cyc  .AND.  bc_ns_cyc ) )  THEN
 
 #if defined( __parallel )   
