Changeset 1314 for palm


Ignore:
Timestamp:
Mar 14, 2014 6:25:17 PM (8 years ago)
Author:
suehring
Message:

Vertical logarithmic interpolation of horizontal particle speed for particles
between roughness height and first vertical grid level.

Location:
palm/trunk/SOURCE
Files:
3 edited

Legend:

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

    r1310 r1314  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Vertical logarithmic interpolation of horizontal particle speed for particles
     23! between roughness height and first vertical grid level.
    2324!
    2425! Former revisions:
     
    6566             e_mean_int, fs_int, lagr_timescale, random_gauss, vv_int
    6667
     68    REAL ::    height_int, height_p, log_z_z0_int, us_int, z_p, d_z_p_z0
     69
    6770    REAL ::  location(1:30,1:3)
    6871    REAL, DIMENSION(1:30) ::  de_dxi, de_dyi, de_dzi, dissi, d_gp_pl, ei
    6972
     73!
     74!-- Determine height of Prandtl layer and distance between Prandtl-layer
     75!-- height and horizontal mean roughness height, which are required for
     76!-- vertical logarithmic interpolation of horizontal particle speeds
     77!-- (for particles below first vertical grid level).
     78    z_p      = zu(nzb+1) - zw(nzb)
     79    d_z_p_z0 = 1.0 / ( z_p - z0_av_global )
    7080
    7181    DO  n = 1, number_of_particles
     
    7585!--    reached
    7686       IF ( ( dt_3d - particles(n)%dt_sum ) < 1E-8 )  CYCLE
    77 
    78 !
    79 !--    Interpolate u velocity-component, determine left, front, bottom
    80 !--    index of u-array
    81        i = ( particles(n)%x + 0.5 * dx ) * ddx
    82        j =   particles(n)%y * ddy
     87!
     88!--    Determine bottom index
    8389       k = ( particles(n)%z + 0.5 * dz * atmos_ocean_sign ) / dz  &
    84            + offset_ocean_nzt                     ! only exact if equidistant
    85 
    86 !
    87 !--    Interpolation of the velocity components in the xy-plane
    88        x  = particles(n)%x + ( 0.5 - i ) * dx
    89        y  = particles(n)%y - j * dy
    90        aa = x**2          + y**2
    91        bb = ( dx - x )**2 + y**2
    92        cc = x**2          + ( dy - y )**2
    93        dd = ( dx - x )**2 + ( dy - y )**2
    94        gg = aa + bb + cc + dd
    95 
    96        u_int_l = ( ( gg - aa ) * u(k,j,i)   + ( gg - bb ) * u(k,j,i+1)   &
    97                  + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * u(k,j+1,i+1) &
    98                  ) / ( 3.0 * gg ) - u_gtrans
    99        IF ( k+1 == nzt+1 )  THEN
    100           u_int = u_int_l
     90              + offset_ocean_nzt       ! only exact if equidistant
     91!
     92!--    Interpolation of the u velocity component onto particle position. 
     93!--    Particles are interpolation bi-linearly in the horizontal and a
     94!--    linearly in the vertical. An exception is made for particles below
     95!--    the first vertical grid level in case of a prandtl layer. In this
     96!--    case the horizontal particle velocity components are determined using
     97!--    Monin-Obukhov relations (if branch).
     98!--    First, check if particle is located below first vertical grid level
     99!--    (Prandtl-layer height)
     100       IF ( prandtl_layer .AND. particles(n)%z < z_p )  THEN
     101!
     102!--       Resolved-scale horizontal particle velocity is zero below z0.
     103          IF ( particles(n)%z < z0_av_global )  THEN
     104
     105             u_int = 0.0
     106
     107          ELSE
     108!
     109!--          Determine the sublayer. Further used as index.
     110             height_p     = ( particles(n)%z - z0_av_global )           &
     111                                  * REAL( number_of_sublayers )         &
     112                                  * d_z_p_z0
     113!
     114!--          Calculate LOG(z/z0) for exact particle height. Therefore,   
     115!--          interpolate linearly between precalculated logarithm.
     116             log_z_z0_int = log_z_z0(INT(height_p))                     &
     117                          + ( height_p - INT(height_p) )                &
     118                          * ( log_z_z0(INT(height_p)+1)                 &
     119                               - log_z_z0(INT(height_p))                &
     120                            )
     121!
     122!--          Neutral solution is applied for all situations, e.g. also for
     123!--          unstable and stable situations. Even though this is not exact
     124!--          this saves a lot of CPU time since several calls of intrinsic
     125!--          FORTRAN procedures (LOG, ATAN) are avoided, This is justified
     126!--          as sensitivity studies revealed no significant effect of
     127!--          using the neutral solution also for un/stable situations.
     128!--          Calculated left and bottom index on u grid.
     129             i      = ( particles(n)%x + 0.5 * dx ) * ddx
     130             j      =   particles(n)%y              * ddy
     131
     132             us_int = 0.5 * ( us(j,i) + us(j,i-1) ) 
     133
     134             u_int  = -usws(j,i) / ( us_int * kappa + 1E-10 )           &
     135                      * log_z_z0_int
     136
     137          ENDIF
     138!
     139!--    Particle above the first grid level. Bi-linear interpolation in the
     140!--    horizontal and linear interpolation in the vertical direction.
    101141       ELSE
    102           u_int_u = ( ( gg-aa ) * u(k+1,j,i)   + ( gg-bb ) * u(k+1,j,i+1)    &
    103                     + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * u(k+1,j+1,i+1)  &
     142!
     143!--       Interpolate u velocity-component, determine left, front, bottom
     144!--       index of u-array. Adopt k index from above
     145          i = ( particles(n)%x + 0.5 * dx ) * ddx
     146          j =   particles(n)%y * ddy
     147!
     148!--       Interpolation of the velocity components in the xy-plane
     149          x  = particles(n)%x + ( 0.5 - i ) * dx
     150          y  = particles(n)%y - j * dy
     151          aa = x**2          + y**2
     152          bb = ( dx - x )**2 + y**2
     153          cc = x**2          + ( dy - y )**2
     154          dd = ( dx - x )**2 + ( dy - y )**2
     155          gg = aa + bb + cc + dd
     156
     157          u_int_l = ( ( gg - aa ) * u(k,j,i)   + ( gg - bb ) * u(k,j,i+1)  &
     158                    + ( gg - cc ) * u(k,j+1,i) + ( gg - dd ) * u(k,j+1,i+1)&
    104159                    ) / ( 3.0 * gg ) - u_gtrans
    105           u_int = u_int_l + ( particles(n)%z - zu(k) ) / dz * &
     160
     161          IF ( k+1 == nzt+1 )  THEN
     162
     163             u_int = u_int_l
     164
     165          ELSE
     166
     167             u_int_u = ( ( gg-aa ) * u(k+1,j,i)   + ( gg-bb ) * u(k+1,j,i+1)   &
     168                       + ( gg-cc ) * u(k+1,j+1,i) + ( gg-dd ) * u(k+1,j+1,i+1) &
     169                       ) / ( 3.0 * gg ) - u_gtrans
     170
     171             u_int = u_int_l + ( particles(n)%z - zu(k) ) / dz *            &
    106172                            ( u_int_u - u_int_l )
     173
     174          ENDIF
     175
    107176       ENDIF
    108177
    109178!
    110 !--    Same procedure for interpolation of the v velocity-component (adopt
    111 !--    index k from u velocity-component)
    112        i =   particles(n)%x * ddx
    113        j = ( particles(n)%y + 0.5 * dy ) * ddy
    114 
    115        x  = particles(n)%x - i * dx
    116        y  = particles(n)%y + ( 0.5 - j ) * dy
    117        aa = x**2          + y**2
    118        bb = ( dx - x )**2 + y**2
    119        cc = x**2          + ( dy - y )**2
    120        dd = ( dx - x )**2 + ( dy - y )**2
    121        gg = aa + bb + cc + dd
    122 
    123        v_int_l = ( ( gg - aa ) * v(k,j,i)   + ( gg - bb ) * v(k,j,i+1)   &
    124                  + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1) &
    125                  ) / ( 3.0 * gg ) - v_gtrans
    126        IF ( k+1 == nzt+1 )  THEN
    127           v_int = v_int_l
     179!--    Same procedure for interpolation of the v velocity-component.
     180       IF ( prandtl_layer .AND. particles(n)%z < z_p )  THEN
     181!
     182!--       Resolved-scale horizontal particle velocity is zero below z0.
     183          IF ( particles(n)%z < z0_av_global )  THEN
     184
     185             v_int = 0.0
     186
     187          ELSE       
     188!
     189!--          Neutral solution is applied for all situations, e.g. also for
     190!--          unstable and stable situations. Even though this is not exact
     191!--          this saves a lot of CPU time since several calls of intrinsic
     192!--          FORTRAN procedures (LOG, ATAN) are avoided, This is justified
     193!--          as sensitivity studies revealed no significant effect of
     194!--          using the neutral solution also for un/stable situations.
     195!--          Calculated left and bottom index on v grid.
     196              i      =   particles(n)%x              * ddx
     197              j      = ( particles(n)%y + 0.5 * dy ) * ddy
     198
     199              us_int = 0.5 * ( us(j,i) + us(j-1,i) )   
     200
     201              v_int  = -vsws(j,i) / ( us_int * kappa + 1E-10 )             &
     202                       * log_z_z0_int
     203
     204          ENDIF
     205!
     206!--    Particle above the first grid level. Bi-linear interpolation in the
     207!--    horizontal and linear interpolation in the vertical direction.
    128208       ELSE
    129           v_int_u = ( ( gg-aa ) * v(k+1,j,i)   + ( gg-bb ) * v(k+1,j,i+1)    &
    130                     + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1)  &
     209          i =   particles(n)%x * ddx
     210          j = ( particles(n)%y + 0.5 * dy ) * ddy
     211          x  = particles(n)%x - i * dx
     212          y  = particles(n)%y + ( 0.5 - j ) * dy
     213          aa = x**2          + y**2
     214          bb = ( dx - x )**2 + y**2
     215          cc = x**2          + ( dy - y )**2
     216          dd = ( dx - x )**2 + ( dy - y )**2
     217          gg = aa + bb + cc + dd
     218
     219          v_int_l = ( ( gg - aa ) * v(k,j,i)   + ( gg - bb ) * v(k,j,i+1)  &
     220                    + ( gg - cc ) * v(k,j+1,i) + ( gg - dd ) * v(k,j+1,i+1)&
    131221                    ) / ( 3.0 * gg ) - v_gtrans
    132           v_int = v_int_l + ( particles(n)%z - zu(k) ) / dz * &
     222          IF ( k+1 == nzt+1 )  THEN
     223             v_int = v_int_l
     224          ELSE
     225             v_int_u = ( ( gg-aa ) * v(k+1,j,i)   + ( gg-bb ) * v(k+1,j,i+1)   &
     226                       + ( gg-cc ) * v(k+1,j+1,i) + ( gg-dd ) * v(k+1,j+1,i+1) &
     227                       ) / ( 3.0 * gg ) - v_gtrans
     228             v_int = v_int_l + ( particles(n)%z - zu(k) ) / dz *           &
    133229                            ( v_int_u - v_int_l )
     230          ENDIF
     231
    134232       ENDIF
    135233
    136234!
    137 !--    Same procedure for interpolation of the w velocity-component (adopt
    138 !--    index i from v velocity-component)
     235!--    Same procedure for interpolation of the w velocity-component
    139236       IF ( vertical_particle_advection(particles(n)%group) )  THEN
     237          i = particles(n)%x * ddx
    140238          j = particles(n)%y * ddy
    141239          k = particles(n)%z / dz + offset_ocean_nzt_m1
  • palm/trunk/SOURCE/lpm_init.f90

    r1310 r1314  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Vertical logarithmic interpolation of horizontal particle speed for particles
     23! between roughness height and first vertical grid level.
    2324!
    2425! Former revisions:
     
    110111    IMPLICIT NONE
    111112
    112     INTEGER ::  i, j, n, nn
     113    INTEGER ::  i, j, k, n, nn
    113114#if defined( __parallel )
    114115    INTEGER, DIMENSION(3) ::  blocklengths, displacements, types
    115116#endif
    116117    LOGICAL ::  uniform_particles_l
    117     REAL    ::  pos_x, pos_y, pos_z
     118    REAL    ::  height_int, height_p, pos_x, pos_y, pos_z, z_p,            &
     119                z0_av_local = 0.0
    118120
    119121
     
    184186                 de_dy(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
    185187                 de_dz(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     188    ENDIF
     189
     190!
     191!-- Allocate array required for logarithmic vertical interpolation of
     192!-- horizontal particle velocities between the surface and the first vertical
     193!-- grid level. In order to avoid repeated CPU cost-intensive CALLS of
     194!-- intrinsic FORTRAN procedure LOG(z/z0), LOG(z/z0) is precalculated for
     195!-- several heights. Splitting into 20 sublayers turned out to be sufficient.
     196!-- To obtain exact height levels of particles, linear interpolation is applied
     197!-- (see lpm_advec.f90).
     198    IF ( prandtl_layer )  THEN
     199       
     200       ALLOCATE ( log_z_z0(0:number_of_sublayers) )
     201       z_p         = zu(nzb+1) - zw(nzb)
     202
     203!
     204!--    Calculate horizontal mean value of z0 used for logartihmic
     205!--    interpolation. Note: this is not exact for heterogeneous z0.
     206!--    However, sensitivity studies showed that the effect is
     207!--    negligible.
     208       z0_av_local  = SUM( z0(nys:nyn,nxl:nxr) )
     209       z0_av_global = 0.0
     210
     211       CALL MPI_ALLREDUCE(z0_av_local, z0_av_global, 1, MPI_REAL, MPI_SUM, &
     212                          comm2d, ierr )
     213
     214       z0_av_global = z0_av_global  / ( ( ny + 1 ) * ( nx + 1 ) )
     215!
     216!--    Horizontal wind speed is zero below and at z0
     217       log_z_z0(0) = 0.0   
     218!
     219!--    Calculate vertical depth of the sublayers
     220       height_int  = ( z_p - z0_av_global ) / REAL( number_of_sublayers )
     221!
     222!--    Precalculate LOG(z/z0)
     223       height_p    = 0.0
     224       DO  k = 1, number_of_sublayers
     225
     226          height_p    = height_p + height_int
     227          log_z_z0(k) = LOG( height_p / z0_av_global )
     228
     229       ENDDO
     230
     231
    186232    ENDIF
    187233
  • palm/trunk/SOURCE/modules.f90

    r1310 r1314  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! + log_z_z0, number_of_sublayers, z0_av_global
    2323!
    2424! Former revisions:
     
    14281428                maximum_number_of_tailpoints = 100,                            &
    14291429                maximum_number_of_tails = 0,                                   &
     1430                number_of_sublayers = 20,                                      &
    14301431                number_of_initial_particles = 0, number_of_particles = 0,      &
    14311432                number_of_particle_groups = 1, number_of_tails = 0,            &
     
    14651466                sgs_wfv_part = 0.3333333, sgs_wfw_part = 0.3333333,            &
    14661467                time_prel = 0.0, time_sort_particles = 0.0,                    &
    1467                 time_write_particle_data = 0.0
     1468                time_write_particle_data = 0.0, z0_av_global
    14681469
    14691470    REAL, DIMENSION(max_number_of_particle_groups) ::  &
     
    14721473                psn = 9999999.9, psr = 9999999.9, pss = 9999999.9,           &
    14731474                pst = 9999999.9, radius = 9999999.9
     1475
     1476    REAL, DIMENSION(:), ALLOCATABLE     ::  log_z_z0
    14741477
    14751478    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  particle_tail_coordinates
Note: See TracChangeset for help on using the changeset viewer.