source: palm/trunk/SOURCE/set_particle_attributes.f90 @ 847

 SUBROUTINE set_particle_attributes

!------------------------------------------------------------------------------!
! Current revisions:
! -----------------
! 
!
! Former revisions:
! -----------------
! $Id: set_particle_attributes.f90 829 2012-02-21 12:24:53Z raasch $
!
! 828 2012-02-21 12:00:36Z raasch
! particle feature color renamed class
!
! 622 2010-12-10 08:08:13Z raasch
! optional barriers included in order to speed up collective operations
!
! 271 2009-03-26 00:47:14Z raasch
! Initial version
!
! Description:
! ------------
! This routine sets certain particle attributes depending on the values that
! other PALM variables have at the current particle position.
!------------------------------------------------------------------------------!

    USE arrays_3d
    USE control_parameters
    USE cpulog
    USE dvrp_variables
    USE grid_variables
    USE indices
    USE interfaces
    USE particle_attributes
    USE pegrid
    USE statistics

    IMPLICIT NONE

    INTEGER ::  i, j, k, n
    REAL    ::  aa, absuv, bb, cc, dd, gg, height, pt_int, pt_int_l, pt_int_u, &
                u_int, u_int_l, u_int_u, v_int, v_int_l, v_int_u, w_int,       &
                w_int_l, w_int_u, x, y


    CALL cpu_log( log_point_s(49), 'set_particle_attrib', 'start' )

!
!-- Set particle color
    IF ( particle_color == 'absuv' )  THEN

!
!--    Set particle color depending on the absolute value of the horizontal
!--    velocity
       DO  n = 1, number_of_particles
!
!--       Interpolate u velocity-component, determine left, front, bottom
!--       index of u-array
          i = ( particles(n)%x + 0.5 * dx ) * ddx
          j =   particles(n)%y * ddy
          k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz  &
              + offset_ocean_nzt                     ! only exact if equidistant

!
!--       Interpolation of the velocity components in the xy-plane
          x  = particles(n)%x + ( 0.5 - i ) * dx
          y  = particles(n)%y - j * dy
          aa = x**2          + y**2
          bb = ( dx - x )**2 + y**2
          cc = x**2          + ( dy - y )**2
          dd = ( dx - x )**2 + ( dy - y )**2
          gg = aa + bb + cc + dd

          u_int_l = ( ( gg - aa ) * u(k,j,i)   + ( gg - bb ) * u(k,j,i+1)   &
                    + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * u(k,j+1,i+1) &
                    ) / ( 3.0 * gg ) - u_gtrans
          IF ( k+1 == nzt+1 )  THEN
             u_int = u_int_l
          ELSE
             u_int_u = ( ( gg-aa ) * u(k+1,j,i)   + ( gg-bb ) * u(k+1,j,i+1)   &
                       + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * u(k+1,j+1,i+1) &
                       ) / ( 3.0 * gg ) - u_gtrans
             u_int = u_int_l + ( particles(n)%z - zu(k) ) / dz * &
                               ( u_int_u - u_int_l )
          ENDIF

!
!--       Same procedure for interpolation of the v velocity-component (adopt
!--       index k from u velocity-component)
          i =   particles(n)%x * ddx
          j = ( particles(n)%y + 0.5 * dy ) * ddy

          x  = particles(n)%x - i * dx
          y  = particles(n)%y + ( 0.5 - j ) * dy
          aa = x**2          + y**2
          bb = ( dx - x )**2 + y**2
          cc = x**2          + ( dy - y )**2
          dd = ( dx - x )**2 + ( dy - y )**2
          gg = aa + bb + cc + dd

          v_int_l = ( ( gg - aa ) * v(k,j,i)   + ( gg - bb ) * v(k,j,i+1)   &
                 + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
                 ) / ( 3.0 * gg ) - v_gtrans
          IF ( k+1 == nzt+1 )  THEN
             v_int = v_int_l
          ELSE
             v_int_u = ( ( gg-aa ) * v(k+1,j,i)   + ( gg-bb ) * v(k+1,j,i+1)   &
                       + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
                       ) / ( 3.0 * gg ) - v_gtrans
             v_int = v_int_l + ( particles(n)%z - zu(k) ) / dz * &
                               ( v_int_u - v_int_l )
          ENDIF

          absuv = SQRT( u_int**2 + v_int**2 )

!
!--       Limit values by the given interval and normalize to interval [0,1]
          absuv = MIN( absuv, color_interval(2) )
          absuv = MAX( absuv, color_interval(1) )

          absuv = ( absuv - color_interval(1) ) / &
                  ( color_interval(2) - color_interval(1) )

!
!--       Number of available colors is defined in init_dvrp
          particles(n)%class = 1 + absuv * ( dvrp_colortable_entries_prt - 1 )

       ENDDO

    ELSEIF ( particle_color == 'pt*' )  THEN
!
!--    Set particle color depending on the resolved scale temperature
!--    fluctuation.
!--    First, calculate the horizontal average of the potential temperature
!--    (This is also done in flow_statistics, but flow_statistics is called
!--    after this routine.)
       sums_l(:,4,0) = 0.0
       DO  i = nxl, nxr
          DO  j =  nys, nyn
             DO  k = nzb, nzt+1
                sums_l(k,4,0) = sums_l(k,4,0) + pt(k,j,i)
             ENDDO
          ENDDO
       ENDDO

#if defined( __parallel )

       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, &
                           MPI_REAL, MPI_SUM, comm2d, ierr )

#else

       sums(:,4) = sums_l(:,4,0)

#endif

       sums(:,4) = sums(:,4) / ngp_2dh(0)

       DO  n = 1, number_of_particles
!
!--       Interpolate temperature to the current particle position
          i = particles(n)%x * ddx
          j = particles(n)%y * ddy
          k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz  &
              + offset_ocean_nzt                     ! only exact if equidistant

          x  = particles(n)%x - i * dx
          y  = particles(n)%y - j * dy
          aa = x**2          + y**2
          bb = ( dx - x )**2 + y**2
          cc = x**2          + ( dy - y )**2
          dd = ( dx - x )**2 + ( dy - y )**2
          gg = aa + bb + cc + dd

          pt_int_l = ( ( gg - aa ) * pt(k,j,i)   + ( gg - bb ) * pt(k,j,i+1)   &
                     + ( gg - cc ) * pt(k,j+1,i) + ( gg - dd ) * pt(k,j+1,i+1) &
                     ) / ( 3.0 * gg ) - sums(k,4)

          pt_int_u = ( ( gg-aa ) * pt(k+1,j,i)   + ( gg-bb ) * pt(k+1,j,i+1)   &
                     + ( gg-cc ) * pt(k+1,j+1,i) + ( gg-dd ) * pt(k+1,j+1,i+1) &
                     ) / ( 3.0 * gg ) - sums(k,4)

          pt_int = pt_int_l + ( particles(n)%z - zu(k) ) / dz * &
                              ( pt_int_u - pt_int_l )

!
!--       Limit values by the given interval and normalize to interval [0,1]
          pt_int = MIN( pt_int, color_interval(2) )
          pt_int = MAX( pt_int, color_interval(1) )

          pt_int = ( pt_int - color_interval(1) ) / &
                   ( color_interval(2) - color_interval(1) )

!
!--       Number of available colors is defined in init_dvrp
          particles(n)%class = 1 + pt_int * ( dvrp_colortable_entries_prt - 1 )

       ENDDO

    ELSEIF ( particle_color == 'z' )  THEN
!
!--    Set particle color depending on the height above the bottom
!--    boundary (z=0)
       DO  n = 1, number_of_particles

          height = particles(n)%z
!
!--       Limit values by the given interval and normalize to interval [0,1]
          height = MIN( height, color_interval(2) )
          height = MAX( height, color_interval(1) )

          height = ( height - color_interval(1) ) / &
                   ( color_interval(2) - color_interval(1) )

!
!--       Number of available colors is defined in init_dvrp
          particles(n)%class = 1 + height * ( dvrp_colortable_entries_prt - 1 )

       ENDDO

    ENDIF

!
!-- Set particle size for dvrp graphics
    IF ( particle_dvrpsize == 'absw' )  THEN

       DO  n = 1, number_of_particles
!
!--       Interpolate w-component to the current particle position
          i = particles(n)%x * ddx
          j = particles(n)%y * ddy
          k = particles(n)%z / dz

          x  = particles(n)%x - i * dx
          y  = particles(n)%y - j * dy
          aa = x**2          + y**2
          bb = ( dx - x )**2 + y**2
          cc = x**2          + ( dy - y )**2
          dd = ( dx - x )**2 + ( dy - y )**2
          gg = aa + bb + cc + dd

          w_int_l = ( ( gg - aa ) * w(k,j,i)   + ( gg - bb ) * w(k,j,i+1)   &
                    + ( gg - cc ) * w(k,j+1,i) + ( gg - dd ) * w(k,j+1,i+1) &
                    ) / ( 3.0 * gg )

          IF ( k+1 == nzt+1 )  THEN
             w_int = w_int_l
          ELSE
             w_int_u = ( ( gg-aa ) * w(k+1,j,i)   + ( gg-bb ) * w(k+1,j,i+1)   &
                       + ( gg-cc ) * w(k+1,j+1,i) + ( gg-dd ) * w(k+1,j+1,i+1) &
                       ) / ( 3.0 * gg )
             w_int = w_int_l + ( particles(n)%z - zw(k) ) / dz * &
                               ( w_int_u - w_int_l )
          ENDIF

!
!--       Limit values by the given interval and normalize to interval [0,1]
          w_int = ABS( w_int )
          w_int = MIN( w_int, dvrpsize_interval(2) )
          w_int = MAX( w_int, dvrpsize_interval(1) )

          w_int = ( w_int - dvrpsize_interval(1) ) / &
                  ( dvrpsize_interval(2) - dvrpsize_interval(1) )

          particles(n)%dvrp_psize = ( 0.25 + w_int * 0.6 ) * dx

       ENDDO

    ENDIF

    CALL cpu_log( log_point_s(49), 'set_particle_attrib', 'stop' )


 END SUBROUTINE set_particle_attributes