SUBROUTINE inflow_turbulence

!------------------------------------------------------------------------------!
! Actual revisions:
! -----------------
! 
!
! Former revisions:
! -----------------
! $Id: inflow_turbulence.f90 151 2008-03-07 13:42:18Z letzel $
!
!
! Description:
! ------------
! Imposing turbulence at the respective inflow using the turbulence
! recycling method of Kataoka and Mizuno (2002).
!------------------------------------------------------------------------------!

    USE arrays_3d
    USE control_parameters
    USE cpulog
    USE grid_variables
    USE indices
    USE interfaces
    USE pegrid


    IMPLICIT NONE

    INTEGER ::  i, imax, j, k, 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

    CALL cpu_log( log_point(40), 'inflow_turbulence', 'start' )

!
!-- Carry out horizontal averaging in the recycling plane
    avpr_l = 0.0
    ngp_pr = ( nzt - nzb + 2 ) * 5
    ngp_ifd = ngp_pr * ( nyn - nys + 3 )

!
!-- First, local averaging within the recycling domain
    IF ( recycling_plane >= nxl )  THEN

       imax = MIN( nxr, recycling_plane )

       DO  i = nxl, imax
          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)

             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 )
!
!-- Now, averaging over all PEs
    CALL MPI_ALLREDUCE( avpr_l(nzb,1), avpr(nzb,1), ngp_pr, MPI_REAL, MPI_SUM, &
                        comm2d, ierr )
#else
    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,*) ' '

!
!-- Calculate the disturbances at the recycling plane
    i = recycling_plane

    IF ( i >= nxl  .AND.  i <= nxr )  THEN

       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
       ENDDO

    ENDIF

!
!-- For parallel runs, send the disturbances to the respective inflow PE
#if defined( __parallel )
    IF ( i >= nxl  .AND.  i <= nxr  .AND.  myid /= id_inflow )  THEN

!       print*, '*** sending id = ', myid, '   send to:', id_inflow
       CALL MPI_SEND( inflow_dist(nzb,nys-1,1), ngp_ifd, MPI_REAL, &
                      id_inflow, 1, comm1dx, ierr )

    ELSEIF ( ( i < nxl  .OR.  i > nxr )  .AND.  myid == id_inflow )  THEN

!       print*, '*** receiving id = ', myid
       CALL MPI_RECV( inflow_dist(nzb,nys-1,1), ngp_ifd, MPI_REAL, &
                      MPI_ANY_SOURCE, 1, comm1dx, status, ierr )

    ENDIF
#endif

!
!-- Add the disturbance at the inflow
    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 )

          ENDDO
       ENDDO

    ENDIF

!
!-- Conserve the volume flow at the inflow in order to avoid generation of
!-- waves in the stable layer
!    IF ( conserve_volume_flow  .AND.  inflow_l )  THEN

!       volume_flow(1)   = 0.0
!       volume_flow_l(1) = 0.0

!       i = 0

!       DO  j = nys, nyn
!
!--       Sum up the volume flow through the south/north boundary
!          DO  k = nzb_2d(j,i) + 1, nzt
!             volume_flow_l(1) = volume_flow_l(1) + u(k,j,i) * dzu(k)
!          ENDDO
!       ENDDO

#if defined( __parallel )   
!       CALL MPI_ALLREDUCE( volume_flow_l(1), volume_flow(1), 1, MPI_REAL, &
!                           MPI_SUM, comm1dy, ierr )   
#else
!       volume_flow = volume_flow_l  
#endif
!       volume_flow_offset(1) = ( volume_flow_initial(1) - volume_flow(1) )    &
!                               / volume_flow_area(1)

!       DO  j = nys-1, nyn+1
!          DO  k = nzb_v_inner(j,i) + 1, nzt
!             u(k,j,i) = u(k,j,i) + volume_flow_offset(1)
!          ENDDO
!       ENDDO

!    ENDIF

    CALL cpu_log( log_point(40), 'inflow_turbulence', 'stop' )


 END SUBROUTINE inflow_turbulence