SUBROUTINE set_particle_attributes !------------------------------------------------------------------------------! ! Current revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: set_particle_attributes.f90 484 2010-02-05 07:36:54Z maronga $ ! ! 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)%color = 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 ) 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)%color = 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)%color = 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