Index: palm/trunk/SOURCE/inflow_turbulence.f90
===================================================================
--- palm/trunk/SOURCE/inflow_turbulence.f90	(revision 623)
+++ palm/trunk/SOURCE/inflow_turbulence.f90	(revision 667)
@@ -4,5 +4,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! Using nbgp recycling planes for a better resolution of the turbulent flow 
+! near the inflow.
 !
 ! Former revisions:
@@ -35,64 +36,68 @@
     IMPLICIT NONE
 
-    INTEGER ::  i, imax, j, k, ngp_ifd, ngp_pr
+    INTEGER ::  i, imax, j, k, l, ngp_ifd, ngp_pr
 
     REAL, DIMENSION(1:2) ::  volume_flow_l, volume_flow_offset
-    REAL, DIMENSION(nzb:nzt+1,5) ::  avpr, avpr_l
-    REAL, DIMENSION(nzb:nzt+1,nys-1:nyn+1,5) ::  inflow_dist
+    REAL, DIMENSION(nzb:nzt+1,5,nbgp) ::  avpr, avpr_l
+    REAL, DIMENSION(nzb:nzt+1,nysg:nyng,5,nbgp) ::  inflow_dist
 
     CALL cpu_log( log_point(40), 'inflow_turbulence', 'start' )
 
 !
-!-- Carry out horizontal averaging in the recycling plane
+!-- Carry out spanwise averaging in the recycling plane
     avpr_l = 0.0
-    ngp_pr = ( nzt - nzb + 2 ) * 5
-    ngp_ifd = ngp_pr * ( nyn - nys + 3 )
+    ngp_pr = ( nzt - nzb + 2 ) * 5 * nbgp
+    ngp_ifd = ngp_pr * ( nyn - nys + 1 + 2 * nbgp )
 
 !
 !-- First, local averaging within the recycling domain
-    IF ( recycling_plane >= nxl )  THEN
-
-       imax = MIN( nxr, recycling_plane )
-
-       DO  i = nxl, imax
+
+    i = recycling_plane
+
+#if defined( __parallel )
+    IF ( myidx == id_recycling )  THEN
+       
+       DO  l = 1, nbgp
           DO  j = nys, nyn
-             DO  k = nzb, nzt+1
-
-                avpr_l(k,1) = avpr_l(k,1) + u(k,j,i)
-                avpr_l(k,2) = avpr_l(k,2) + v(k,j,i)
-                avpr_l(k,3) = avpr_l(k,3) + w(k,j,i)
-                avpr_l(k,4) = avpr_l(k,4) + pt(k,j,i)
-                avpr_l(k,5) = avpr_l(k,5) + e(k,j,i)
+             DO  k = nzb, nzt + 1
+
+                avpr_l(k,1,l) = avpr_l(k,1,l) + u(k,j,i)
+                avpr_l(k,2,l) = avpr_l(k,2,l) + v(k,j,i)
+                avpr_l(k,3,l) = avpr_l(k,3,l) + w(k,j,i)
+                avpr_l(k,4,l) = avpr_l(k,4,l) + pt(k,j,i)
+                avpr_l(k,5,l) = avpr_l(k,5,l) + e(k,j,i)
 
              ENDDO
           ENDDO
-       ENDDO
-
-    ENDIF
-
-!    WRITE (9,*) '*** averaged profiles avpr_l'
-!    DO  k = nzb, nzt+1
-!       WRITE (9,'(F5.1,1X,F5.1,1X,F5.1,1X,F6.1,1X,F7.2)') avpr_l(k,1),avpr_l(k,2),avpr_l(k,3),avpr_l(k,4),avpr_l(k,5)
-!    ENDDO
-!    WRITE (9,*) ' '
-
-#if defined( __parallel )
+          i = i + 1
+       ENDDO
+
+    ENDIF
 !
 !-- Now, averaging over all PEs
     IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-    CALL MPI_ALLREDUCE( avpr_l(nzb,1), avpr(nzb,1), ngp_pr, MPI_REAL, MPI_SUM, &
-                        comm2d, ierr )
+    CALL MPI_ALLREDUCE( avpr_l(nzb,1,1), avpr(nzb,1,1), ngp_pr,  & 
+                    MPI_REAL, MPI_SUM, comm2d, ierr )
+
 #else
+    DO  l = 1, nbgp
+       DO  j = nys, nyn
+          DO  k = nzb, nzt + 1
+
+             avpr_l(k,1,l) = avpr_l(k,1,l) + u(k,j,i)
+             avpr_l(k,2,l) = avpr_l(k,2,l) + v(k,j,i)
+             avpr_l(k,3,l) = avpr_l(k,3,l) + w(k,j,i)
+             avpr_l(k,4,l) = avpr_l(k,4,l) + pt(k,j,i)
+             avpr_l(k,5,l) = avpr_l(k,5,l) + e(k,j,i)
+
+          ENDDO
+       ENDDO
+       i = i + 1 
+    ENDDO
+    
     avpr = avpr_l
 #endif
 
-    avpr = avpr / ( ( ny + 1 ) * ( recycling_plane + 1 ) )
-
-!    WRITE (9,*) '*** averaged profiles'
-!    DO  k = nzb, nzt+1
-!       WRITE (9,'(F5.1,1X,F5.1,1X,F5.1,1X,F6.1,1X,F7.2)') avpr(k,1),avpr(k,2),avpr(k,3),avpr(k,4),avpr(k,5)
-!    ENDDO
-!    WRITE (9,*) ' '
-
+    avpr = avpr / ( ny + 1 )
 !
 !-- Calculate the disturbances at the recycling plane
@@ -101,29 +106,34 @@
 #if defined( __parallel )
     IF ( myidx == id_recycling )  THEN
-
-       DO  j = nys-1, nyn+1
+       DO  l = 1, nbgp
+          DO  j = nysg, nyng
+             DO  k = nzb, nzt + 1
+
+                inflow_dist(k,j,1,l) = u(k,j,i+1) - avpr(k,1,l)
+                inflow_dist(k,j,2,l) = v(k,j,i)   - avpr(k,2,l)
+                inflow_dist(k,j,3,l) = w(k,j,i)   - avpr(k,3,l)
+                inflow_dist(k,j,4,l) = pt(k,j,i)  - avpr(k,4,l)
+                inflow_dist(k,j,5,l) = e(k,j,i)   - avpr(k,5,l)
+              
+            ENDDO
+          ENDDO
+          i = i + 1
+       ENDDO
+
+    ENDIF
+#else
+    DO  l = 1, nbgp
+       DO  j = nysg, nyng
           DO  k = nzb, nzt+1
 
-              inflow_dist(k,j,1) = u(k,j,i+1) - avpr(k,1)
-              inflow_dist(k,j,2) = v(k,j,i)   - avpr(k,2)
-              inflow_dist(k,j,3) = w(k,j,i)   - avpr(k,3)
-              inflow_dist(k,j,4) = pt(k,j,i)  - avpr(k,4)
-              inflow_dist(k,j,5) = e(k,j,i)   - avpr(k,5)
-
-          ENDDO
-       ENDDO
-
-    ENDIF
-#else
-    DO  j = nys-1, nyn+1
-       DO  k = nzb, nzt+1
-
-          inflow_dist(k,j,1) = u(k,j,i+1) - avpr(k,1)
-          inflow_dist(k,j,2) = v(k,j,i)   - avpr(k,2)
-          inflow_dist(k,j,3) = w(k,j,i)   - avpr(k,3)
-          inflow_dist(k,j,4) = pt(k,j,i)  - avpr(k,4)
-          inflow_dist(k,j,5) = e(k,j,i)   - avpr(k,5)
-
-       ENDDO
+             inflow_dist(k,j,1,l) = u(k,j,i+1) - avpr(k,1,l)
+             inflow_dist(k,j,2,l) = v(k,j,i)   - avpr(k,2,l)
+             inflow_dist(k,j,3,l) = w(k,j,i)   - avpr(k,3,l)
+             inflow_dist(k,j,4,l) = pt(k,j,i)  - avpr(k,4,l)
+             inflow_dist(k,j,5,l) = e(k,j,i)   - avpr(k,5,l)
+              
+          ENDDO
+       ENDDO
+       i = i + 1
     ENDDO
 #endif
@@ -134,5 +144,5 @@
     IF ( myidx == id_recycling  .AND.  myidx /= id_inflow )  THEN
 
-       CALL MPI_SEND( inflow_dist(nzb,nys-1,1), ngp_ifd, MPI_REAL, &
+       CALL MPI_SEND( inflow_dist(nzb,nysg,1,1), ngp_ifd, MPI_REAL, &
                       id_inflow, 1, comm1dx, ierr )
 
@@ -140,5 +150,5 @@
 
        inflow_dist = 0.0
-       CALL MPI_RECV( inflow_dist(nzb,nys-1,1), ngp_ifd, MPI_REAL, &
+       CALL MPI_RECV( inflow_dist(nzb,nysg,1,1), ngp_ifd, MPI_REAL, &
                       id_recycling, 1, comm1dx, status, ierr )
 
@@ -150,28 +160,17 @@
     IF ( nxl == 0 )  THEN
 
-       DO  j = nys-1, nyn+1
-          DO  k = nzb, nzt+1
-
-!              WRITE (9,*) 'j=',j,' k=',k
-!              WRITE (9,*) 'mean_u = ', mean_inflow_profiles(k,1), ' dist_u = ',&
-!                          inflow_dist(k,j,1)
-!              WRITE (9,*) 'mean_v = ', mean_inflow_profiles(k,2), ' dist_v = ',&
-!                          inflow_dist(k,j,2)
-!              WRITE (9,*) 'mean_w = 0.0', ' dist_w = ',&
-!                          inflow_dist(k,j,3)
-!              WRITE (9,*) 'mean_pt = ', mean_inflow_profiles(k,4), ' dist_pt = ',&
-!                          inflow_dist(k,j,4)
-!              WRITE (9,*) 'mean_e = ', mean_inflow_profiles(k,5), ' dist_e = ',&
-!                          inflow_dist(k,j,5)
-              u(k,j,0)   = mean_inflow_profiles(k,1) + &
-                           inflow_dist(k,j,1) * inflow_damping_factor(k)
-              v(k,j,-1)  = mean_inflow_profiles(k,2) + &
-                           inflow_dist(k,j,2) * inflow_damping_factor(k)
-              w(k,j,-1)  = inflow_dist(k,j,3) * inflow_damping_factor(k)
-              pt(k,j,-1) = mean_inflow_profiles(k,4) + &
-                           inflow_dist(k,j,4) * inflow_damping_factor(k)
-              e(k,j,-1)  = mean_inflow_profiles(k,5) + &
-                           inflow_dist(k,j,5) * inflow_damping_factor(k)
-              e(k,j,-1)  = MAX( e(k,j,-1), 0.0 )
+       DO  j = nysg, nyng
+          DO  k = nzb, nzt + 1
+
+              u(k,j,-nbgp+1:0)   = mean_inflow_profiles(k,1) + &
+                           inflow_dist(k,j,1,1:nbgp) * inflow_damping_factor(k)
+              v(k,j,-nbgp:-1)  = mean_inflow_profiles(k,2) + &
+                           inflow_dist(k,j,2,1:nbgp) * inflow_damping_factor(k)
+              w(k,j,-nbgp:-1)  = inflow_dist(k,j,3,1:nbgp) * inflow_damping_factor(k)
+              pt(k,j,-nbgp:-1) = mean_inflow_profiles(k,4) + &
+                           inflow_dist(k,j,4,1:nbgp) * inflow_damping_factor(k)
+              e(k,j,-nbgp:-1)  = mean_inflow_profiles(k,5) + &
+                           inflow_dist(k,j,5,1:nbgp) * inflow_damping_factor(k)
+              e(k,j,-nbgp:-1)  = MAX( e(k,j,-nbgp:-1), 0.0 )
 
           ENDDO
