Ignore:
Timestamp:
Mar 25, 2009 7:43:59 AM (15 years ago)
Author:
raasch
Message:

further updates concerning new dvr-features + documentation

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/set_particle_attributes.f90

    r264 r266  
    3232
    3333    INTEGER ::  i, j, k, n
    34     REAL    ::  aa, bb, cc, dd, gg, pt_int, pt_int_l, pt_int_u, w_int, &
    35                 w_int_l, w_int_u, x, y
     34    REAL    ::  aa, absuv, bb, cc, dd, gg, pt_int, pt_int_l, pt_int_u, u_int, &
     35                u_int_l, u_int_u, v_int, v_int_l, v_int_u, w_int, w_int_l,    &
     36                w_int_u, x, y
    3637
    3738
     
    4142!-- Set particle color
    4243    IF ( particle_color == 'absuv' )  THEN
     44
     45!
     46!--    Set particle color depending on the absolute value of the horizontal
     47!--    velocity
     48       DO  n = 1, number_of_particles
     49!
     50!--       Interpolate u velocity-component, determine left, front, bottom
     51!--       index of u-array
     52          i = ( particles(n)%x + 0.5 * dx ) * ddx
     53          j =   particles(n)%y * ddy
     54          k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz  &
     55              + offset_ocean_nzt                     ! only exact if equidistant
     56
     57!
     58!--       Interpolation of the velocity components in the xy-plane
     59          x  = particles(n)%x + ( 0.5 - i ) * dx
     60          y  = particles(n)%y - j * dy
     61          aa = x**2          + y**2
     62          bb = ( dx - x )**2 + y**2
     63          cc = x**2          + ( dy - y )**2
     64          dd = ( dx - x )**2 + ( dy - y )**2
     65          gg = aa + bb + cc + dd
     66
     67          u_int_l = ( ( gg - aa ) * u(k,j,i)   + ( gg - bb ) * u(k,j,i+1)   &
     68                    + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * u(k,j+1,i+1) &
     69                    ) / ( 3.0 * gg ) - u_gtrans
     70          IF ( k+1 == nzt+1 )  THEN
     71             u_int = u_int_l
     72          ELSE
     73             u_int_u = ( ( gg-aa ) * u(k+1,j,i)   + ( gg-bb ) * u(k+1,j,i+1)   &
     74                       + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * u(k+1,j+1,i+1) &
     75                       ) / ( 3.0 * gg ) - u_gtrans
     76             u_int = u_int_l + ( particles(n)%z - zu(k) ) / dz * &
     77                               ( u_int_u - u_int_l )
     78          ENDIF
     79
     80!
     81!--       Same procedure for interpolation of the v velocity-component (adopt
     82!--       index k from u velocity-component)
     83          i =   particles(n)%x * ddx
     84          j = ( particles(n)%y + 0.5 * dy ) * ddy
     85
     86          x  = particles(n)%x - i * dx
     87          y  = particles(n)%y + ( 0.5 - j ) * dy
     88          aa = x**2          + y**2
     89          bb = ( dx - x )**2 + y**2
     90          cc = x**2          + ( dy - y )**2
     91          dd = ( dx - x )**2 + ( dy - y )**2
     92          gg = aa + bb + cc + dd
     93
     94          v_int_l = ( ( gg - aa ) * v(k,j,i)   + ( gg - bb ) * v(k,j,i+1)   &
     95                 + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
     96                 ) / ( 3.0 * gg ) - v_gtrans
     97          IF ( k+1 == nzt+1 )  THEN
     98             v_int = v_int_l
     99          ELSE
     100             v_int_u = ( ( gg-aa ) * v(k+1,j,i)   + ( gg-bb ) * v(k+1,j,i+1)   &
     101                       + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
     102                       ) / ( 3.0 * gg ) - v_gtrans
     103             v_int = v_int_l + ( particles(n)%z - zu(k) ) / dz * &
     104                               ( v_int_u - v_int_l )
     105          ENDIF
     106
     107          absuv = SQRT( u_int**2 + v_int**2 )
     108
     109!
     110!--       Limit values by the given interval and normalize to interval [0,1]
     111          absuv = MIN( absuv, color_interval(2) )
     112          absuv = MAX( absuv, color_interval(1) )
     113
     114          absuv = ( absuv - color_interval(1) ) / &
     115                  ( color_interval(2) - color_interval(1) )
     116
     117!
     118!--       Number of available colors is defined in init_dvrp
     119          particles(n)%color = 1 + absuv * ( dvrp_colortable_entries_prt - 1 )
     120
     121       ENDDO
    43122
    44123    ELSEIF ( particle_color == 'pt*' )  THEN
     
    110189          particles(n)%color = 1 + pt_int * ( dvrp_colortable_entries_prt - 1 )
    111190
    112     ENDDO
     191       ENDDO
    113192
    114193    ELSEIF ( particle_color == 'z' )  THEN
     194!
     195!--    Set particle color depending on the height above the bottom
     196!--    boundary (z=0)
     197       DO  n = 1, number_of_particles
     198
     199          height = particles(n)%z
     200!
     201!--       Limit values by the given interval and normalize to interval [0,1]
     202          height = MIN( height, color_interval(2) )
     203          height = MAX( height, color_interval(1) )
     204
     205          height = ( height - color_interval(1) ) / &
     206                   ( color_interval(2) - color_interval(1) )
     207
     208!
     209!--       Number of available colors is defined in init_dvrp
     210          particles(n)%color = 1 + height * ( dvrp_colortable_entries_prt - 1 )
     211
     212       ENDDO
    115213
    116214    ENDIF
     
    149247          ENDIF
    150248
    151           particles(n)%dvrp_psize = ( 0.2 + MIN( 0.6, ABS( w_int ) / 6.0 ) ) * &
    152                                     dx * 1.4
     249!
     250!--       Limit values by the given interval and normalize to interval [0,1]
     251          w_int = ABS( w_int )
     252          w_int = MIN( w_int, dvrpsize_interval(2) )
     253          w_int = MAX( w_int, dvrpsize_interval(1) )
     254
     255          w_int = ( w_int - dvrpsize_interval(1) ) / &
     256                  ( dvrpsize_interval(2) - dvrpsize_interval(1) )
     257
     258          particles(n)%dvrp_psize = ( 0.25 + w_int * 0.6 ) * dx
    153259
    154260       ENDDO
Note: See TracChangeset for help on using the changeset viewer.