Ignore:
Timestamp:
Dec 14, 2017 6:46:24 PM (4 years ago)
Author:
suehring
Message:

Particle reflections at downward-facing walls; revision of particle speed interpolations at walls; bugfixes in get_topography_index and in date constants

File:
1 edited

Legend:

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

    r2696 r2698  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Particle reflections at downward-facing walls implemented. Moreover,
     23! reflections are adjusted to revised particle grid box location.
     24! (responsible Philipp Thiele)
    2325!
    2426! Former revisions:
     
    9294!> (see offset_ocean_*)
    9395!------------------------------------------------------------------------------!
    94  SUBROUTINE lpm_boundary_conds( location )
     96 SUBROUTINE lpm_boundary_conds( location , i, j, k )
    9597 
    9698
     
    108110
    109111    USE indices,                                                               &
    110         ONLY:  nxl, nxr, nyn, nys, nz, nzb
     112        ONLY:  nxl, nxr, nyn, nys, nz, nzb, wall_flags_0,nyng,nysg
    111113
    112114    USE kinds
     
    119121    USE pegrid
    120122
    121     USE surface_mod,                                                           &
    122         ONLY:  get_topography_top_index
    123 
    124123    IMPLICIT NONE
    125124
    126125    CHARACTER (LEN=*) ::  location     !<
     126   
     127    INTEGER(iwp), INTENT(IN) ::  i !<
     128    INTEGER(iwp), INTENT(IN) ::  j !<
     129    INTEGER(iwp), INTENT(IN) ::  k !<
    127130   
    128131    INTEGER(iwp) ::  inc            !< dummy for sorting algorithmus
     
    133136    INTEGER(iwp) ::  jr             !< dummy for sorting algorithmus
    134137    INTEGER(iwp) ::  j1             !< grid index (y) of old particle position
    135     INTEGER(iwp) ::  j2             !< grid index (x) of current particle position
    136     INTEGER(iwp) ::  j3             !< grid index (x) of intermediate particle position
    137     INTEGER(iwp) ::  k_wall         !< vertical index of topography top
     138    INTEGER(iwp) ::  j2             !< grid index (y) of current particle position
     139    INTEGER(iwp) ::  j3             !< grid index (y) of intermediate particle position
     140    INTEGER(iwp) ::  k1             !< grid index (z) of old particle position
     141    INTEGER(iwp) ::  k2             !< grid index (z) of current particle position
     142    INTEGER(iwp) ::  k3             !< grid index (z) of intermediate particle position
    138143    INTEGER(iwp) ::  n              !< particle number
    139144    INTEGER(iwp) ::  t_index        !< running index for intermediate particle timesteps in reflection algorithmus
    140145    INTEGER(iwp) ::  t_index_number !< number of intermediate particle timesteps in reflection algorithmus
    141     INTEGER(iwp) ::  tmp_x          !< dummy for sorting algorithmus
    142     INTEGER(iwp) ::  tmp_y          !< dummy for sorting algorithmus
     146    INTEGER(iwp) ::  tmp_x          !< dummy for sorting algorithm
     147    INTEGER(iwp) ::  tmp_y          !< dummy for sorting algorithm
     148    INTEGER(iwp) ::  tmp_z          !< dummy for sorting algorithm
    143149
    144150    INTEGER(iwp), DIMENSION(0:10) :: x_ind(0:10) = 0 !< index array (x) of intermediate particle positions
    145     INTEGER(iwp), DIMENSION(0:10) :: y_ind(0:10) = 0 !< index array (x) of intermediate particle positions
     151    INTEGER(iwp), DIMENSION(0:10) :: y_ind(0:10) = 0 !< index array (y) of intermediate particle positions
     152    INTEGER(iwp), DIMENSION(0:10) :: z_ind(0:10) = 0 !< index array (z) of intermediate particle positions
    146153   
    147154    LOGICAL  ::  cross_wall_x    !< flag to check if particle reflection along x is necessary
    148155    LOGICAL  ::  cross_wall_y    !< flag to check if particle reflection along y is necessary
    149     LOGICAL  ::  downwards       !< flag to check if particle reflection along z is necessary (only if particle move downwards)
     156    LOGICAL  ::  cross_wall_z    !< flag to check if particle reflection along z is necessary
    150157    LOGICAL  ::  reflect_x       !< flag to check if particle is already reflected along x
    151158    LOGICAL  ::  reflect_y       !< flag to check if particle is already reflected along y
     
    156163    LOGICAL  ::  x_wall_reached  !< flag to check if particle has already reached wall
    157164    LOGICAL  ::  y_wall_reached  !< flag to check if particle has already reached wall
     165    LOGICAL  ::  z_wall_reached  !< flag to check if particle has already reached wall
    158166
    159167    LOGICAL, DIMENSION(0:10) ::  reach_x  !< flag to check if particle is at a yz-wall
     
    177185    REAL(wp) ::  xwall          !< location of wall in x
    178186    REAL(wp) ::  ywall          !< location of wall in y
    179     REAL(wp) ::  zwall1         !< location of wall in z (old grid box)
    180     REAL(wp) ::  zwall2         !< location of wall in z (current grid box)
    181     REAL(wp) ::  zwall3         !< location of wall in z (old y, current x)
    182     REAL(wp) ::  zwall4         !< location of wall in z (current y, old x)
     187    REAL(wp) ::  zwall          !< location of wall in z
    183188
    184189    REAL(wp), DIMENSION(0:10) ::  t  !< reflection time
     
    263268          i2 = particles(n)%x * ddx
    264269          j2 = particles(n)%y * ddy
     270          IF (zw(k)   < particles(n)%z ) k2 = k + 1
     271          IF (zw(k-1) > particles(n)%z ) k2 = k - 1
    265272!
    266273!--       Save current particle positions
     
    275282!
    276283!--       Obtain x/y indices for old particle positions
    277           i1 = pos_x_old * ddx
    278           j1 = pos_y_old * ddy
     284          i1 = i
     285          j1 = j
     286          k1 = k
    279287!
    280288!--       Determine horizontal as well as vertical walls at which particle can
     
    283291!--       Wall to the right
    284292          IF ( prt_x > pos_x_old )  THEN
    285              xwall = ( i1 + 0.5_wp ) * dx
     293             xwall = ( i1 + 1 ) * dx
    286294!
    287295!--       Wall to the left
    288296          ELSE
    289              xwall = ( i1 - 0.5_wp ) * dx
     297             xwall = i1 * dx
    290298          ENDIF
    291299!
     
    293301!--       Wall to the north
    294302          IF ( prt_y > pos_y_old )  THEN
    295              ywall = ( j1 + 0.5_wp ) * dy
     303             ywall = ( j1 +1 ) * dy
    296304!--       Wall to the south
    297305          ELSE
    298              ywall = ( j1 - 0.5_wp ) * dy
    299           ENDIF
    300 !
    301 !--       Walls aligned in xy layer at which particle can be possiblly reflected.
    302 !--       The construct of MERGE and BTEST is used to determine the topography-
    303 !--       top index (former nzb_s_inner).
    304           zwall1 = zw( get_topography_top_index( j2, i2, 's' ) )                                             
    305           zwall2 = zw( get_topography_top_index( j1, i1, 's' ) )
    306           zwall3 = zw( get_topography_top_index( j1, i2, 's' ) )
    307           zwall4 = zw( get_topography_top_index( j2, i1, 's' ) )
     306             ywall = j1 * dy
     307          ENDIF
     308
     309          IF ( prt_z > pos_z_old ) THEN
     310             zwall = zw(k)
     311          ELSE
     312             zwall = zw(k-1)
     313          ENDIF     
    308314!
    309315!--       Initialize flags to check if particle reflection is necessary
    310           downwards    = .FALSE.
    311316          cross_wall_x = .FALSE.
    312317          cross_wall_y = .FALSE.
     318          cross_wall_z = .FALSE.
    313319!
    314320!--       Initialize flags to check if a wall is reached
     
    318324!
    319325!--       Initialize flags to check if a particle was already reflected
    320           reflect_x = .FALSE.
    321           reflect_y = .FALSE.
    322           reflect_z = .FALSE.
    323 !
    324 !--       Initialize flags to check if a vertical wall is already crossed.
     326          reflect_x    = .FALSE.
     327          reflect_y    = .FALSE.
     328          reflect_z    = .FALSE.
     329!
     330!--       Initialize flags to check if a wall is already crossed.
    325331!--       ( Required to obtain correct indices. )
    326332          x_wall_reached = .FALSE.
    327333          y_wall_reached = .FALSE.
     334          z_wall_reached = .FALSE.
    328335!
    329336!--       Initialize time array
     
    342349          x_ind(t_index)   = i2
    343350          y_ind(t_index)   = j1
     351          z_ind(t_index)   = k1
    344352          reach_x(t_index) = .TRUE.
    345353          reach_y(t_index) = .FALSE.
     
    360368          x_ind(t_index)   = i1
    361369          y_ind(t_index)   = j2
     370          z_ind(t_index)   = k1
    362371          reach_x(t_index) = .FALSE.
    363372          reach_y(t_index) = .TRUE.
     
    369378!
    370379!--       z-direction
    371 !--       At first, check if particle moves downwards. Only in this case a
    372 !--       particle can be reflected vertically.
    373           IF ( prt_z < pos_z_old )  THEN
    374 
    375              downwards = .TRUE.
    376              dum       =  1.0_wp / MERGE( MAX( prt_z - pos_z_old,  1E-30_wp ), &
    377                                           MIN( prt_z - pos_z_old, -1E-30_wp ), &
    378                                           prt_z > pos_z_old )
    379 
    380              t(t_index)       = ( zwall1 - pos_z_old ) * dum
    381              x_ind(t_index)   = i2
    382              y_ind(t_index)   = j2
    383              reach_x(t_index) = .FALSE.
    384              reach_y(t_index) = .FALSE.
    385              reach_z(t_index) = .TRUE.
    386              IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )            &
    387                 t_index = t_index + 1
    388 
    389              reach_x(t_index) = .FALSE.
    390              reach_y(t_index) = .FALSE.
    391              reach_z(t_index) = .TRUE.
    392              t(t_index)       = ( zwall2 - pos_z_old ) * dum
    393              x_ind(t_index)   = i1
    394              y_ind(t_index)   = j1
    395              IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )            &
    396                 t_index = t_index + 1
    397 
    398              reach_x(t_index) = .FALSE.
    399              reach_y(t_index) = .FALSE.
    400              reach_z(t_index) = .TRUE.
    401              t(t_index)       = ( zwall3 - pos_z_old ) * dum
    402              x_ind(t_index)   = i2
    403              y_ind(t_index)   = j1
    404              IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )            &
    405                 t_index = t_index + 1
    406 
    407              reach_x(t_index) = .FALSE.
    408              reach_y(t_index) = .FALSE.
    409              reach_z(t_index) = .TRUE.
    410              t(t_index)       = ( zwall4 - pos_z_old ) * dum
    411              x_ind(t_index)   = i1
    412              y_ind(t_index)   = j2
    413              IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )            &
    414                 t_index = t_index + 1
    415 
    416           END IF
     380          t(t_index) = (zwall - pos_z_old )                                    &
     381                     / MERGE( MAX( prt_z - pos_z_old,  1E-30_wp ),             &
     382                              MIN( prt_z - pos_z_old, -1E-30_wp ),             &
     383                              prt_z > pos_z_old )
     384                     
     385          x_ind(t_index)   = i1
     386          y_ind(t_index)   = j1
     387          z_ind(t_index)   = k2
     388          reach_x(t_index) = .FALSE.
     389          reach_y(t_index) = .FALSE.
     390          reach_z(t_index) = .TRUE.
     391          IF( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp) THEN
     392             t_index      = t_index + 1
     393             cross_wall_z = .TRUE.
     394          ENDIF
     395         
    417396          t_index_number = t_index - 1
    418397!
    419398!--       Carry out reflection only if particle reaches any wall
    420           IF ( cross_wall_x .OR. cross_wall_y .OR. downwards )  THEN
     399          IF ( cross_wall_x .OR. cross_wall_y .OR. cross_wall_z )  THEN
    421400!
    422401!--          Sort fractional timesteps in ascending order. Also sort the
     
    435414                   tmp_x       = x_ind(ir)
    436415                   tmp_y       = y_ind(ir)
     416                   tmp_z       = z_ind(ir)
    437417                   tmp_reach_x = reach_x(ir)
    438418                   tmp_reach_y = reach_y(ir)
     
    443423                      x_ind(jr)   = x_ind(jr-inc)
    444424                      y_ind(jr)   = y_ind(jr-inc)
     425                      z_ind(jr)   = z_ind(jr-inc)
    445426                      reach_x(jr) = reach_x(jr-inc)
    446427                      reach_y(jr) = reach_y(jr-inc)
     
    452433                   x_ind(jr)   = tmp_x
    453434                   y_ind(jr)   = tmp_y
     435                   z_ind(jr)   = tmp_z
    454436                   reach_x(jr) = tmp_reach_x
    455437                   reach_y(jr) = tmp_reach_y
     
    480462                i3 = x_ind(t_index)
    481463                j3 = y_ind(t_index)
     464                k3 = z_ind(t_index)
    482465!
    483466!--             Check which wall is already reached
    484467                IF ( .NOT. x_wall_reached )  x_wall_reached = reach_x(t_index)
    485                 IF ( .NOT. y_wall_reached )  y_wall_reached = reach_y(t_index)
     468                IF ( .NOT. y_wall_reached )  y_wall_reached = reach_y(t_index)
     469                IF ( .NOT. z_wall_reached )  z_wall_reached = reach_z(t_index)
    486470!
    487471!--             Check if a particle needs to be reflected at any yz-wall. If
    488472!--             necessary, carry out reflection. Please note, a security
    489 !--             constant is required, as the particle position do not
     473!--             constant is required, as the particle position does not
    490474!--             necessarily exactly match the wall location due to rounding
    491 !--             errors. At first, determine index of topography top at (j3,i3) 
    492                 k_wall = get_topography_top_index( j3, i3, 's' )
    493                 IF ( ABS( pos_x - xwall ) < eps      .AND.                     &
    494                      pos_z <= zw(k_wall)             .AND.                     &
    495                      reach_x(t_index)                .AND.                     &
     475!--             errors.
     476                IF ( reach_x(t_index)                      .AND.               &
     477                     ABS( pos_x - xwall ) < eps            .AND.               &
     478                     .NOT. BTEST(wall_flags_0(k3,j3,i3),0) .AND.               &
    496479                     .NOT. reflect_x )  THEN
    497 !
     480!
     481!
    498482!--                Reflection in x-direction.
    499483!--                Ensure correct reflection by MIN/MAX functions, depending on
    500484!--                direction of particle transport.
    501 !--                Due to rounding errors pos_x do not exactly matches the wall
     485!--                Due to rounding errors pos_x does not exactly match the wall
    502486!--                location, leading to erroneous reflection.             
    503487                   pos_x = MERGE( MIN( 2.0_wp * xwall - pos_x, xwall ),        &
     
    508492                   particles(n)%speed_x = - particles(n)%speed_x
    509493!
    510 !--                Change also sign of subgrid-scale particle speed
     494!--                Also change sign of subgrid-scale particle speed
    511495                   particles(n)%rvar1 = - particles(n)%rvar1
    512496!
     
    514498                   reflect_x          = .TRUE.
    515499!
    516 !--                As particle do not crosses any further yz-wall during
     500!--                As the particle does not cross any further yz-wall during
    517501!--                this timestep, set further x-indices to the current one.
    518502                   x_ind(t_index:t_index_number) = i1
     
    522506                ELSEIF ( x_wall_reached .AND. .NOT. reflect_x )  THEN
    523507                    x_ind(t_index:t_index_number) = i2
    524                 ENDIF
     508                ENDIF !particle reflection in x direction done
     509
    525510!
    526511!--             Check if a particle needs to be reflected at any xz-wall. If
    527 !--             necessary, carry out reflection. At first, determine index of
    528 !--             topography top at (j3,i3) 
    529                 k_wall = get_topography_top_index( j3, i3, 's' )
    530                 IF ( ABS( pos_y - ywall ) < eps      .AND.                     &
    531                      pos_z <= zw(k_wall)             .AND.                     &
    532                      reach_y(t_index)                .AND.                     &
    533                      .NOT. reflect_y ) THEN
    534 
     512!--             necessary, carry out reflection. Please note, a security
     513!--             constant is required, as the particle position does not
     514!--             necessarily exactly match the wall location due to rounding
     515!--             errors.
     516                WRITE(9,*) 'wall_test_y nysg:',nysg ,' y: ',particles(n)%y,' j3: ',j3,'nyng:', nyng 
     517                FLUSH(9)
     518                WRITE(9,*) BTEST(wall_flags_0(k3,j3,i3),0)
     519                FLUSH(9)
     520                IF ( reach_y(t_index)                      .AND.               &
     521                     ABS( pos_y - ywall ) < eps            .AND.               &
     522                     .NOT. BTEST(wall_flags_0(k3,j3,i3),0) .AND.               &
     523                     .NOT. reflect_y )  THEN
     524!
     525!
     526!--                Reflection in y-direction.
     527!--                Ensure correct reflection by MIN/MAX functions, depending on
     528!--                direction of particle transport.
     529!--                Due to rounding errors pos_y does not exactly match the wall
     530!--                location, leading to erroneous reflection.             
    535531                   pos_y = MERGE( MIN( 2.0_wp * ywall - pos_y, ywall ),        &
    536532                                  MAX( 2.0_wp * ywall - pos_y, ywall ),        &
    537                                   particles(n)%y > ywall )
    538 
    539                    particles(n)%speed_y = - particles(n)%speed_y     
    540                    particles(n)%rvar2   = - particles(n)%rvar2       
    541      
    542                    reflect_y            = .TRUE.
     533                                  particles(n)%y > ywall )
     534!
     535!--                Change sign of particle speed                     
     536                   particles(n)%speed_y = - particles(n)%speed_y
     537!
     538!--                Also change sign of subgrid-scale particle speed
     539                   particles(n)%rvar2 = - particles(n)%rvar2
     540!
     541!--                Set flag that reflection along y is already done
     542                   reflect_y          = .TRUE.
     543!
     544!--                As the particle does not cross any further xz-wall during
     545!--                this timestep, set further y-indices to the current one.
    543546                   y_ind(t_index:t_index_number) = j1
    544 
     547!
     548!--             If particle already reached the wall but was not reflected,
     549!--             set further y-indices to the new one.
    545550                ELSEIF ( y_wall_reached .AND. .NOT. reflect_y )  THEN
    546                    y_ind(t_index:t_index_number) = j2
    547                 ENDIF
     551                    y_ind(t_index:t_index_number) = j2
     552                ENDIF !particle reflection in y direction done
     553               
    548554!
    549555!--             Check if a particle needs to be reflected at any xy-wall. If
    550 !--             necessary, carry out reflection.
    551                 IF ( downwards .AND. reach_z(t_index) .AND.                    &
     556!--             necessary, carry out reflection. Please note, a security
     557!--             constant is required, as the particle position does not
     558!--             necessarily exactly match the wall location due to rounding
     559!--             errors.
     560                IF ( reach_z(t_index)                      .AND.               &
     561                     ABS( pos_z - zwall ) < eps            .AND.               &
     562                     .NOT. BTEST(wall_flags_0(k3,j3,i3),0) .AND.               &
    552563                     .NOT. reflect_z )  THEN
    553 !
    554 !--                Determine index of topography top at (j3,i3) and chick if
    555 !--                particle is below. 
    556                    k_wall = get_topography_top_index( j3, i3, 's' )
    557                    IF ( pos_z - zw(k_wall) < eps ) THEN
    558  
    559                       pos_z = MAX( 2.0_wp * zw(k_wall) - pos_z,    &
    560                                    zw(k_wall) )
    561 
    562                       particles(n)%speed_z = - particles(n)%speed_z
    563                       particles(n)%rvar3   = - particles(n)%rvar3
    564 
    565                       reflect_z            = .TRUE.
    566 
    567                    ENDIF
    568 
    569                 ENDIF
     564!
     565!
     566!--                Reflection in z-direction.
     567!--                Ensure correct reflection by MIN/MAX functions, depending on
     568!--                direction of particle transport.
     569!--                Due to rounding errors pos_z does not exactly match the wall
     570!--                location, leading to erroneous reflection.             
     571                   pos_z = MERGE( MIN( 2.0_wp * zwall - pos_z, zwall ),        &
     572                                  MAX( 2.0_wp * zwall - pos_z, zwall ),        &
     573                                  particles(n)%z > zwall )
     574!
     575!--                Change sign of particle speed                     
     576                   particles(n)%speed_z = - particles(n)%speed_z
     577!
     578!--                Also change sign of subgrid-scale particle speed
     579                   particles(n)%rvar3 = - particles(n)%rvar3
     580!
     581!--                Set flag that reflection along z is already done
     582                   reflect_z          = .TRUE.
     583!
     584!--                As the particle does not cross any further xy-wall during
     585!--                this timestep, set further z-indices to the current one.
     586                   z_ind(t_index:t_index_number) = k1
     587!
     588!--             If particle already reached the wall but was not reflected,
     589!--             set further z-indices to the new one.
     590                ELSEIF ( z_wall_reached .AND. .NOT. reflect_z )  THEN
     591                    z_ind(t_index:t_index_number) = k2
     592                ENDIF !particle reflection in z direction done               
     593               
    570594!
    571595!--             Swap time
Note: See TracChangeset for help on using the changeset viewer.