Ignore:
Timestamp:
Jun 9, 2016 4:25:25 PM (5 years ago)
Author:
suehring
Message:

several bugfixes in particle model and serial mode

File:
1 edited

Legend:

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

    r1823 r1929  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Rewritten wall reflection
    2222!
    2323! Former revisions:
     
    102102    CHARACTER (LEN=*) ::  range     !<
    103103   
    104     INTEGER(iwp) ::  i              !<
    105     INTEGER(iwp) ::  inc            !<
    106     INTEGER(iwp) ::  ir             !<
    107     INTEGER(iwp) ::  i1             !<
    108     INTEGER(iwp) ::  i2             !<
    109     INTEGER(iwp) ::  i3             !<
    110     INTEGER(iwp) ::  i5             !<
    111     INTEGER(iwp) ::  j              !<
    112     INTEGER(iwp) ::  jr             !<
    113     INTEGER(iwp) ::  j1             !<
    114     INTEGER(iwp) ::  j2             !<
    115     INTEGER(iwp) ::  j3             !<
    116     INTEGER(iwp) ::  j5             !<
    117     INTEGER(iwp) ::  k              !<
    118     INTEGER(iwp) ::  k1             !<
    119     INTEGER(iwp) ::  k2             !<
    120     INTEGER(iwp) ::  k3             !<
    121     INTEGER(iwp) ::  k5             !<
    122     INTEGER(iwp) ::  n              !<
    123     INTEGER(iwp) ::  t_index        !<
    124     INTEGER(iwp) ::  t_index_number !<
     104    INTEGER(iwp) ::  inc            !< dummy for sorting algorithmus
     105    INTEGER(iwp) ::  ir             !< dummy for sorting algorithmus
     106    INTEGER(iwp) ::  i1             !< grid index (x) of old particle position
     107    INTEGER(iwp) ::  i2             !< grid index (x) of current particle position
     108    INTEGER(iwp) ::  i3             !< grid index (x) of intermediate particle position
     109    INTEGER(iwp) ::  jr             !< dummy for sorting algorithmus
     110    INTEGER(iwp) ::  j1             !< grid index (y) of old particle position
     111    INTEGER(iwp) ::  j2             !< grid index (x) of current particle position
     112    INTEGER(iwp) ::  j3             !< grid index (x) of intermediate particle position
     113    INTEGER(iwp) ::  n              !< particle number
     114    INTEGER(iwp) ::  t_index        !< running index for intermediate particle timesteps in reflection algorithmus
     115    INTEGER(iwp) ::  t_index_number !< number of intermediate particle timesteps in reflection algorithmus
     116    INTEGER(iwp) ::  tmp_x          !< dummy for sorting algorithmus
     117    INTEGER(iwp) ::  tmp_y          !< dummy for sorting algorithmus
     118
     119    INTEGER(iwp), DIMENSION(0:10) :: x_ind(0:10) = 0 !< index array (x) of intermediate particle positions
     120    INTEGER(iwp), DIMENSION(0:10) :: y_ind(0:10) = 0 !< index array (x) of intermediate particle positions
    125121   
    126     LOGICAL  ::  reflect_x   !<
    127     LOGICAL  ::  reflect_y   !<
    128     LOGICAL  ::  reflect_z   !<
    129 
    130     REAL(wp) ::  dt_particle !<
    131     REAL(wp) ::  pos_x       !<
    132     REAL(wp) ::  pos_x_old   !<
    133     REAL(wp) ::  pos_y       !<
    134     REAL(wp) ::  pos_y_old   !<
    135     REAL(wp) ::  pos_z       !<
    136     REAL(wp) ::  pos_z_old   !<
    137     REAL(wp) ::  prt_x       !<
    138     REAL(wp) ::  prt_y       !<
    139     REAL(wp) ::  prt_z       !<
    140     REAL(wp) ::  t(1:200)    !<
    141     REAL(wp) ::  tmp_t       !<
    142     REAL(wp) ::  xline       !<
    143     REAL(wp) ::  yline       !<
    144     REAL(wp) ::  zline       !<
     122    LOGICAL  ::  cross_wall_x    !< flag to check if particle reflection along x is necessary
     123    LOGICAL  ::  cross_wall_y    !< flag to check if particle reflection along y is necessary
     124    LOGICAL  ::  downwards       !< flag to check if particle reflection along z is necessary (only if particle move downwards)
     125    LOGICAL  ::  reflect_x       !< flag to check if particle is already reflected along x
     126    LOGICAL  ::  reflect_y       !< flag to check if particle is already reflected along y
     127    LOGICAL  ::  reflect_z       !< flag to check if particle is already reflected along z
     128    LOGICAL  ::  tmp_reach_x     !< dummy for sorting algorithmus
     129    LOGICAL  ::  tmp_reach_y     !< dummy for sorting algorithmus
     130    LOGICAL  ::  tmp_reach_z     !< dummy for sorting algorithmus
     131    LOGICAL  ::  x_wall_reached  !< flag to check if particle has already reached wall
     132    LOGICAL  ::  y_wall_reached  !< flag to check if particle has already reached wall
     133
     134    LOGICAL, DIMENSION(0:10) ::  reach_x  !< flag to check if particle is at a yz-wall
     135    LOGICAL, DIMENSION(0:10) ::  reach_y  !< flag to check if particle is at a xz-wall
     136    LOGICAL, DIMENSION(0:10) ::  reach_z  !< flag to check if particle is at a xy-wall
     137
     138    REAL(wp) ::  dt_particle    !< particle timestep
     139    REAL(wp) ::  dum            !< dummy argument
     140    REAL(wp) ::  eps = 1E-10_wp !< security number to check if particle has reached a wall
     141    REAL(wp) ::  pos_x          !< intermediate particle position (x)
     142    REAL(wp) ::  pos_x_old      !< particle position (x) at previous particle timestep
     143    REAL(wp) ::  pos_y          !< intermediate particle position (y)
     144    REAL(wp) ::  pos_y_old      !< particle position (y) at previous particle timestep
     145    REAL(wp) ::  pos_z          !< intermediate particle position (z)
     146    REAL(wp) ::  pos_z_old      !< particle position (z) at previous particle timestep
     147    REAL(wp) ::  prt_x          !< current particle position (x)
     148    REAL(wp) ::  prt_y          !< current particle position (y)
     149    REAL(wp) ::  prt_z          !< current particle position (z)
     150    REAL(wp) ::  t_old          !< previous reflection time
     151    REAL(wp) ::  tmp_t          !< dummy for sorting algorithmus
     152    REAL(wp) ::  xwall          !< location of wall in x
     153    REAL(wp) ::  ywall          !< location of wall in y
     154    REAL(wp) ::  zwall1         !< location of wall in z (old grid box)
     155    REAL(wp) ::  zwall2         !< location of wall in z (current grid box)
     156    REAL(wp) ::  zwall3         !< location of wall in z (old y, current x)
     157    REAL(wp) ::  zwall4         !< location of wall in z (current y, old x)
     158
     159    REAL(wp), DIMENSION(0:10) ::  t  !< reflection time
     160
    145161
    146162    IF ( range == 'bottom/top' )  THEN
     
    211227    ELSEIF ( range == 'walls' )  THEN
    212228
     229
    213230       CALL cpu_log( log_point_s(48), 'lpm_wall_reflect', 'start' )
    214231
    215        reflect_x = .FALSE.
    216        reflect_y = .FALSE.
    217        reflect_z = .FALSE.
    218 
    219232       DO  n = 1, number_of_particles
    220 
     233!
     234!--       Recalculate particle timestep
    221235          dt_particle = particles(n)%age - particles(n)%age_m
    222 
     236!
     237!--       Obtain x/y indices for current particle position
    223238          i2 = ( particles(n)%x + 0.5_wp * dx ) * ddx
    224239          j2 = ( particles(n)%y + 0.5_wp * dy ) * ddy
    225           k2 = particles(n)%z / dz + 1 + offset_ocean_nzt_m1
    226 
     240!
     241!--       Save current particle positions
    227242          prt_x = particles(n)%x
    228243          prt_y = particles(n)%y
    229244          prt_z = particles(n)%z
    230 
    231 !
    232 !--       If the particle position is below the surface, it has to be reflected
    233           IF ( k2 <= nzb_s_inner(j2,i2)  .AND.  nzb_s_inner(j2,i2) /=0 )  THEN
    234 
    235              pos_x_old = particles(n)%x - particles(n)%speed_x * dt_particle
    236              pos_y_old = particles(n)%y - particles(n)%speed_y * dt_particle
    237              pos_z_old = particles(n)%z - particles(n)%speed_z * dt_particle
    238              i1 = ( pos_x_old + 0.5_wp * dx ) * ddx
    239              j1 = ( pos_y_old + 0.5_wp * dy ) * ddy
    240              k1 = pos_z_old / dz + offset_ocean_nzt_m1
    241 
    242 !
    243 !--          Case 1
    244              IF ( particles(n)%x > pos_x_old  .AND.  particles(n)%y > pos_y_old )&
    245              THEN
    246                 t_index = 1
    247 
    248                 DO  i = i1, i2
    249                    xline      = i * dx + 0.5_wp * dx
    250                    t(t_index) = ( xline - pos_x_old ) / &
    251                                 ( particles(n)%x - pos_x_old )
    252                    t_index    = t_index + 1
     245!
     246!--       Recalculate old particle positions
     247          pos_x_old = particles(n)%x - particles(n)%speed_x * dt_particle
     248          pos_y_old = particles(n)%y - particles(n)%speed_y * dt_particle
     249          pos_z_old = particles(n)%z - particles(n)%speed_z * dt_particle
     250!
     251!--       Obtain x/y indices for old particle positions
     252          i1 = ( pos_x_old + 0.5_wp * dx ) * ddx
     253          j1 = ( pos_y_old + 0.5_wp * dy ) * ddy
     254!
     255!--       Determine horizontal as well as vertical walls at which particle can
     256!--       be potentially reflected.
     257!--       Start with walls aligned in yz layer.
     258!--       Wall to the right
     259          IF ( prt_x > pos_x_old )  THEN
     260             xwall = ( i1 + 0.5_wp ) * dx
     261!
     262!--       Wall to the left
     263          ELSE
     264             xwall = ( i1 - 0.5_wp ) * dx
     265          ENDIF
     266!
     267!--       Walls aligned in xz layer
     268!--       Wall to the north
     269          IF ( prt_y > pos_y_old )  THEN
     270             ywall = ( j1 + 0.5_wp ) * dy
     271!--       Wall to the south
     272          ELSE
     273             ywall = ( j1 - 0.5_wp ) * dy
     274          ENDIF
     275!
     276!--       Walls aligned in xy layer at which particle can be possiblly reflected
     277          zwall1 = zw(nzb_s_inner(j2,i2))
     278          zwall2 = zw(nzb_s_inner(j1,i1))
     279          zwall3 = zw(nzb_s_inner(j1,i2))
     280          zwall4 = zw(nzb_s_inner(j2,i1))
     281!
     282!--       Initialize flags to check if particle reflection is necessary
     283          downwards    = .FALSE.
     284          cross_wall_x = .FALSE.
     285          cross_wall_y = .FALSE.
     286!
     287!--       Initialize flags to check if a wall is reached
     288          reach_x      = .FALSE.
     289          reach_y      = .FALSE.
     290          reach_z      = .FALSE.
     291!
     292!--       Initialize flags to check if a particle was already reflected
     293          reflect_x = .FALSE.
     294          reflect_y = .FALSE.
     295          reflect_z = .FALSE.
     296!
     297!--       Initialize flags to check if a vertical wall is already crossed.
     298!--       ( Required to obtain correct indices. )
     299          x_wall_reached = .FALSE.
     300          y_wall_reached = .FALSE.
     301!
     302!--       Initialize time array
     303          t     = 0.0_wp
     304!
     305!--       Check if particle can reach any wall. This case, calculate the
     306!--       fractional time needed to reach this wall. Store this fractional
     307!--       timestep in array t. Moreover, store indices for these grid
     308!--       boxes where the respective wall belongs to. 
     309!--       Start with x-direction.
     310          t_index    = 1
     311          t(t_index) = ( xwall - pos_x_old )                                   &
     312                     / MERGE( MAX( prt_x - pos_x_old,  1E-30_wp ),             &
     313                              MIN( prt_x - pos_x_old, -1E-30_wp ),             &
     314                              prt_x > pos_x_old )
     315          x_ind(t_index)   = i2
     316          y_ind(t_index)   = j1
     317          reach_x(t_index) = .TRUE.
     318          reach_y(t_index) = .FALSE.
     319          reach_z(t_index) = .FALSE.
     320!
     321!--       Store these values only if particle really reaches any wall. t must
     322!--       be in a interval between [0:1].
     323          IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )  THEN
     324             t_index      = t_index + 1
     325             cross_wall_x = .TRUE.
     326          ENDIF
     327!
     328!--       y-direction
     329          t(t_index) = ( ywall - pos_y_old )                                   &
     330                     / MERGE( MAX( prt_y - pos_y_old,  1E-30_wp ),             &
     331                              MIN( prt_y - pos_y_old, -1E-30_wp ),             &
     332                              prt_y > pos_y_old )
     333          x_ind(t_index)   = i1
     334          y_ind(t_index)   = j2
     335          reach_x(t_index) = .FALSE.
     336          reach_y(t_index) = .TRUE.
     337          reach_z(t_index) = .FALSE.
     338          IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )  THEN
     339             t_index      = t_index + 1
     340             cross_wall_y = .TRUE.
     341          ENDIF
     342!
     343!--       z-direction
     344!--       At first, check if particle moves downwards. Only in this case a
     345!--       particle can be reflected vertically.
     346          IF ( prt_z < pos_z_old )  THEN
     347
     348             downwards = .TRUE.
     349             dum       =  1.0_wp / MERGE( MAX( prt_z - pos_z_old,  1E-30_wp ), &
     350                                          MIN( prt_z - pos_z_old, -1E-30_wp ), &
     351                                          prt_z > pos_z_old )
     352
     353             t(t_index)       = ( zwall1 - pos_z_old ) * dum
     354             x_ind(t_index)   = i2
     355             y_ind(t_index)   = j2
     356             reach_x(t_index) = .FALSE.
     357             reach_y(t_index) = .FALSE.
     358             reach_z(t_index) = .TRUE.
     359             IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )            &
     360                t_index = t_index + 1
     361
     362             reach_x(t_index) = .FALSE.
     363             reach_y(t_index) = .FALSE.
     364             reach_z(t_index) = .TRUE.
     365             t(t_index)       = ( zwall2 - pos_z_old ) * dum
     366             x_ind(t_index)   = i1
     367             y_ind(t_index)   = j1
     368             IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )            &
     369                t_index = t_index + 1
     370
     371             reach_x(t_index) = .FALSE.
     372             reach_y(t_index) = .FALSE.
     373             reach_z(t_index) = .TRUE.
     374             t(t_index)       = ( zwall3 - pos_z_old ) * dum
     375             x_ind(t_index)   = i2
     376             y_ind(t_index)   = j1
     377             IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )            &
     378                t_index = t_index + 1
     379
     380             reach_x(t_index) = .FALSE.
     381             reach_y(t_index) = .FALSE.
     382             reach_z(t_index) = .TRUE.
     383             t(t_index)       = ( zwall4 - pos_z_old ) * dum
     384             x_ind(t_index)   = i1
     385             y_ind(t_index)   = j2
     386             IF ( t(t_index) <= 1.0_wp .AND. t(t_index) >= 0.0_wp )            &
     387                t_index = t_index + 1
     388
     389          END IF
     390          t_index_number = t_index - 1
     391!
     392!--       Carry out reflection only if particle reaches any wall
     393          IF ( cross_wall_x .OR. cross_wall_y .OR. downwards )  THEN
     394!
     395!--          Sort fractional timesteps in ascending order. Also sort the
     396!--          corresponding indices and flag according to the time interval a 
     397!--          particle reaches the respective wall.
     398             inc = 1
     399             jr  = 1
     400             DO WHILE ( inc <= t_index_number )
     401                inc = 3 * inc + 1
     402             ENDDO
     403
     404             DO WHILE ( inc > 1 )
     405                inc = inc / 3
     406                DO  ir = inc+1, t_index_number
     407                   tmp_t       = t(ir)
     408                   tmp_x       = x_ind(ir)
     409                   tmp_y       = y_ind(ir)
     410                   tmp_reach_x = reach_x(ir)
     411                   tmp_reach_y = reach_y(ir)
     412                   tmp_reach_z = reach_z(ir)
     413                   jr    = ir
     414                   DO WHILE ( t(jr-inc) > tmp_t )
     415                      t(jr)       = t(jr-inc)
     416                      x_ind(jr)   = x_ind(jr-inc)
     417                      y_ind(jr)   = y_ind(jr-inc)
     418                      reach_x(jr) = reach_x(jr-inc)
     419                      reach_y(jr) = reach_y(jr-inc)
     420                      reach_z(jr) = reach_z(jr-inc)
     421                      jr    = jr - inc
     422                      IF ( jr <= inc )  EXIT
     423                   ENDDO
     424                   t(jr)       = tmp_t
     425                   x_ind(jr)   = tmp_x
     426                   y_ind(jr)   = tmp_y
     427                   reach_x(jr) = tmp_reach_x
     428                   reach_y(jr) = tmp_reach_y
     429                   reach_z(jr) = tmp_reach_z
    253430                ENDDO
    254 
    255                 DO  j = j1, j2
    256                    yline      = j * dy + 0.5_wp * dy
    257                    t(t_index) = ( yline - pos_y_old ) / &
    258                                 ( particles(n)%y - pos_y_old )
    259                    t_index    = t_index + 1
    260                 ENDDO
    261 
    262                 IF ( particles(n)%z < pos_z_old )  THEN
    263                    DO  k = k1, k2 , -1
    264                       zline      = k * dz
    265                       t(t_index) = ( pos_z_old - zline ) / &
    266                                    ( pos_z_old - particles(n)%z )
    267                       t_index    = t_index + 1
    268                    ENDDO
     431             ENDDO
     432!
     433!--          Initialize temporary particle positions
     434             pos_x = pos_x_old
     435             pos_y = pos_y_old
     436             pos_z = pos_z_old
     437!
     438!--          Loop over all times a particle possibly moves into a new grid box
     439             t_old = 0.0_wp
     440             DO t_index = 1, t_index_number
     441!           
     442!--             Calculate intermediate particle position according to the
     443!--             timesteps a particle reaches any wall.
     444                pos_x = pos_x + ( t(t_index) - t_old ) * dt_particle           &
     445                                                       * particles(n)%speed_x
     446                pos_y = pos_y + ( t(t_index) - t_old ) * dt_particle           &
     447                                                       * particles(n)%speed_y
     448                pos_z = pos_z + ( t(t_index) - t_old ) * dt_particle           &
     449                                                       * particles(n)%speed_z
     450!
     451!--             Obtain x/y grid indices for intermediate particle position from
     452!--             sorted index array
     453                i3 = x_ind(t_index)
     454                j3 = y_ind(t_index)
     455!
     456!--             Check which wall is already reached
     457                IF ( .NOT. x_wall_reached )  x_wall_reached = reach_x(t_index)
     458                IF ( .NOT. y_wall_reached )  y_wall_reached = reach_y(t_index)
     459!
     460!--             Check if a particle needs to be reflected at any yz-wall. If
     461!--             necessary, carry out reflection. Please note, a security
     462!--             constant is required, as the particle position do not
     463!--             necessarily exactly match the wall location due to rounding
     464!--             errors.   
     465                IF ( ABS( pos_x - xwall ) < eps      .AND.                     &
     466                     pos_z <= zw(nzb_s_inner(j3,i3)) .AND.                     &
     467                     reach_x(t_index)                .AND.                     &
     468                     .NOT. reflect_x )  THEN
     469!
     470!--                Reflection in x-direction.
     471!--                Ensure correct reflection by MIN/MAX functions, depending on
     472!--                direction of particle transport.
     473!--                Due to rounding errors pos_x do not exactly matches the wall
     474!--                location, leading to erroneous reflection.             
     475                   pos_x = MERGE( MIN( 2.0_wp * xwall - pos_x, xwall ),        &
     476                                  MAX( 2.0_wp * xwall - pos_x, xwall ),        &
     477                                  particles(n)%x > xwall )
     478!
     479!--                Change sign of particle speed                     
     480                   particles(n)%speed_x = - particles(n)%speed_x
     481!
     482!--                Change also sign of subgrid-scale particle speed
     483                   particles(n)%rvar1 = - particles(n)%rvar1
     484!
     485!--                Set flag that reflection along x is already done
     486                   reflect_x          = .TRUE.
     487!
     488!--                As particle do not crosses any further yz-wall during
     489!--                this timestep, set further x-indices to the current one.
     490                   x_ind(t_index:t_index_number) = i1
     491!
     492!--             If particle already reached the wall but was not reflected,
     493!--             set further x-indices to the new one.
     494                ELSEIF ( x_wall_reached .AND. .NOT. reflect_x )  THEN
     495                    x_ind(t_index:t_index_number) = i2
    269496                ENDIF
    270 
    271                 t_index_number = t_index - 1
    272 
    273 !
    274 !--             Sorting t
    275                 inc = 1
    276                 jr  = 1
    277                 DO WHILE ( inc <= t_index_number )
    278                    inc = 3 * inc + 1
    279                 ENDDO
    280 
    281                 DO WHILE ( inc > 1 )
    282                    inc = inc / 3
    283                    DO  ir = inc+1, t_index_number
    284                       tmp_t = t(ir)
    285                       jr    = ir
    286                       DO WHILE ( t(jr-inc) > tmp_t )
    287                          t(jr) = t(jr-inc)
    288                          jr    = jr - inc
    289                          IF ( jr <= inc )  EXIT
    290                       ENDDO
    291                       t(jr) = tmp_t
    292                    ENDDO
    293                 ENDDO
    294 
    295          case1: DO  t_index = 1, t_index_number
    296 
    297                    pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
    298                    pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
    299                    pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
    300 
    301                    i3 = ( pos_x + 0.5_wp * dx ) * ddx   
    302                    j3 = ( pos_y + 0.5_wp * dy ) * ddy
    303                    k3 = pos_z / dz + offset_ocean_nzt_m1
    304 
    305                    i5 = pos_x * ddx
    306                    j5 = pos_y * ddy
    307                    k5 = pos_z / dz + offset_ocean_nzt_m1
    308 
    309                    IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
    310                         nzb_s_inner(j5,i3) /= 0 )  THEN
    311 
    312                       IF ( pos_z == nzb_s_inner(j5,i3) * dz  .AND. &
    313                            k3 == nzb_s_inner(j5,i3) )  THEN
    314                          reflect_z = .TRUE.
    315                          EXIT case1
    316                       ENDIF
     497!
     498!--             Check if a particle needs to be reflected at any xz-wall. If
     499!--             necessary, carry out reflection.
     500                IF ( ABS( pos_y - ywall ) < eps      .AND.                     &
     501                     pos_z <= zw(nzb_s_inner(j3,i3)) .AND.                     &
     502                     reach_y(t_index)                .AND.                     &
     503                     .NOT. reflect_y ) THEN
     504
     505                   pos_y = MERGE( MIN( 2.0_wp * ywall - pos_y, ywall ),        &
     506                                  MAX( 2.0_wp * ywall - pos_y, ywall ),        &
     507                                  particles(n)%y > ywall )
     508
     509                   particles(n)%speed_y = - particles(n)%speed_y     
     510                   particles(n)%rvar2   = - particles(n)%rvar2       
     511     
     512                   reflect_y            = .TRUE.
     513                   y_ind(t_index:t_index_number) = j1
     514
     515                ELSEIF ( y_wall_reached .AND. .NOT. reflect_y )  THEN
     516                   y_ind(t_index:t_index_number) = j2
     517                ENDIF
     518!
     519!--             Check if a particle needs to be reflected at any xy-wall. If
     520!--             necessary, carry out reflection.
     521                IF ( downwards .AND. reach_z(t_index) .AND.                    &
     522                     .NOT. reflect_z )  THEN
     523                   IF ( pos_z - zw(nzb_s_inner(j3,i3)) < eps ) THEN
     524 
     525                      pos_z = MAX( 2.0_wp * zw(nzb_s_inner(j3,i3)) - pos_z,    &
     526                                   zw(nzb_s_inner(j3,i3)) )
     527
     528                      particles(n)%speed_z = - particles(n)%speed_z
     529                      particles(n)%rvar3   = - particles(n)%rvar3
     530
     531                      reflect_z            = .TRUE.
    317532
    318533                   ENDIF
    319534
    320                    IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
    321                         nzb_s_inner(j3,i5) /= 0 )  THEN
    322 
    323                       IF ( pos_z == nzb_s_inner(j3,i5) * dz  .AND. &
    324                            k3 == nzb_s_inner(j3,i5) )  THEN
    325                          reflect_z = .TRUE.
    326                          EXIT case1
    327                       ENDIF
    328 
    329                    ENDIF
    330 
    331                    IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
    332                         nzb_s_inner(j3,i3) /= 0 )  THEN
    333 
    334                       IF ( pos_z == nzb_s_inner(j3,i3) * dz  .AND. &
    335                            k3 == nzb_s_inner(j3,i3) )  THEN
    336                          reflect_z = .TRUE.
    337                          EXIT case1
    338                       ENDIF
    339 
    340                       IF ( pos_y == ( j3 * dy - 0.5_wp * dy )  .AND. &
    341                            pos_z < nzb_s_inner(j3,i3) * dz )  THEN
    342                          reflect_y = .TRUE.
    343                          EXIT case1
    344                       ENDIF
    345 
    346                       IF ( pos_x == ( i3 * dx - 0.5_wp * dx )  .AND. &
    347                            pos_z < nzb_s_inner(j3,i3) * dz )  THEN
    348                          reflect_x = .TRUE.
    349                          EXIT case1
    350                       ENDIF
    351 
    352                    ENDIF
    353 
    354                 ENDDO case1
    355 
    356 !
    357 !--          Case 2
    358              ELSEIF ( particles(n)%x > pos_x_old  .AND. &
    359                       particles(n)%y < pos_y_old )  THEN
    360 
    361                 t_index = 1
    362 
    363                 DO  i = i1, i2
    364                    xline      = i * dx + 0.5_wp * dx
    365                    t(t_index) = ( xline - pos_x_old ) / &
    366                                 ( particles(n)%x - pos_x_old )
    367                    t_index    = t_index + 1
    368                 ENDDO
    369 
    370                 DO  j = j1, j2, -1
    371                    yline      = j * dy - 0.5_wp * dy
    372                    t(t_index) = ( pos_y_old - yline ) / &
    373                                 ( pos_y_old - particles(n)%y )
    374                    t_index    = t_index + 1
    375                 ENDDO
    376 
    377                 IF ( particles(n)%z < pos_z_old )  THEN
    378                    DO  k = k1, k2 , -1
    379                       zline      = k * dz
    380                       t(t_index) = ( pos_z_old - zline ) / &
    381                                    ( pos_z_old - particles(n)%z )
    382                       t_index    = t_index + 1
    383                    ENDDO
    384535                ENDIF
    385                 t_index_number = t_index-1
    386 
    387 !
    388 !--             Sorting t
    389                 inc = 1
    390                 jr  = 1
    391                 DO WHILE ( inc <= t_index_number )
    392                    inc = 3 * inc + 1
    393                 ENDDO
    394 
    395                 DO WHILE ( inc > 1 )
    396                    inc = inc / 3
    397                    DO  ir = inc+1, t_index_number
    398                       tmp_t = t(ir)
    399                       jr    = ir
    400                       DO WHILE ( t(jr-inc) > tmp_t )
    401                          t(jr) = t(jr-inc)
    402                          jr    = jr - inc
    403                          IF ( jr <= inc )  EXIT
    404                       ENDDO
    405                       t(jr) = tmp_t
    406                    ENDDO
    407                 ENDDO
    408 
    409          case2: DO  t_index = 1, t_index_number
    410 
    411                    pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
    412                    pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
    413                    pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
    414 
    415                    i3 = ( pos_x + 0.5_wp * dx ) * ddx
    416                    j3 = ( pos_y + 0.5_wp * dy ) * ddy
    417                    k3 = pos_z / dz + offset_ocean_nzt_m1
    418 
    419                    i5 = pos_x * ddx
    420                    j5 = pos_y * ddy
    421                    k5 = pos_z / dz + offset_ocean_nzt_m1
    422 
    423                    IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
    424                         nzb_s_inner(j3,i5) /= 0 )  THEN
    425 
    426                       IF ( pos_z == nzb_s_inner(j3,i5) * dz  .AND. &
    427                            k3 == nzb_s_inner(j3,i5) )  THEN
    428                          reflect_z = .TRUE.
    429                          EXIT case2
    430                       ENDIF
    431 
    432                    ENDIF
    433 
    434                    IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
    435                         nzb_s_inner(j3,i3) /= 0 )  THEN
    436 
    437                       IF ( pos_z == nzb_s_inner(j3,i3) * dz  .AND. &
    438                            k3 == nzb_s_inner(j3,i3) )  THEN
    439                          reflect_z = .TRUE.
    440                          EXIT case2
    441                       ENDIF
    442 
    443                       IF ( pos_x == ( i3 * dx - 0.5_wp * dx )  .AND. &
    444                            pos_z < nzb_s_inner(j3,i3) * dz )  THEN
    445                          reflect_x = .TRUE.
    446                          EXIT case2
    447                       ENDIF
    448 
    449                    ENDIF
    450 
    451                    IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
    452                         nzb_s_inner(j5,i3) /= 0 )  THEN
    453 
    454                       IF ( pos_z == nzb_s_inner(j5,i3) * dz  .AND. &
    455                            k3 == nzb_s_inner(j5,i3) )  THEN
    456                          reflect_z = .TRUE.
    457                          EXIT case2
    458                       ENDIF
    459 
    460                       IF ( pos_y == ( j5 * dy + 0.5_wp * dy )  .AND. &
    461                            pos_z < nzb_s_inner(j5,i3) * dz )  THEN
    462                          reflect_y = .TRUE.
    463                          EXIT case2
    464                       ENDIF
    465 
    466                    ENDIF
    467 
    468                 ENDDO case2
    469 
    470 !
    471 !--          Case 3
    472              ELSEIF ( particles(n)%x < pos_x_old  .AND. &
    473                       particles(n)%y > pos_y_old )  THEN
    474 
    475                 t_index = 1
    476 
    477                 DO  i = i1, i2, -1
    478                    xline      = i * dx - 0.5_wp * dx
    479                    t(t_index) = ( pos_x_old - xline ) / &
    480                                 ( pos_x_old - particles(n)%x )
    481                    t_index    = t_index + 1
    482                 ENDDO
    483 
    484                 DO  j = j1, j2
    485                    yline      = j * dy + 0.5_wp * dy
    486                    t(t_index) = ( yline - pos_y_old ) / &
    487                                 ( particles(n)%y - pos_y_old )
    488                    t_index    = t_index + 1
    489                 ENDDO
    490 
    491                 IF ( particles(n)%z < pos_z_old )  THEN
    492                    DO  k = k1, k2 , -1
    493                       zline      = k * dz
    494                       t(t_index) = ( pos_z_old - zline ) / &
    495                                    ( pos_z_old - particles(n)%z )
    496                       t_index    = t_index + 1
    497                    ENDDO
    498                 ENDIF
    499                 t_index_number = t_index - 1
    500 
    501 !
    502 !--             Sorting t
    503                 inc = 1
    504                 jr  = 1
    505 
    506                 DO WHILE ( inc <= t_index_number )
    507                    inc = 3 * inc + 1
    508                 ENDDO
    509 
    510                 DO WHILE ( inc > 1 )
    511                    inc = inc / 3
    512                    DO  ir = inc+1, t_index_number
    513                       tmp_t = t(ir)
    514                       jr    = ir
    515                       DO WHILE ( t(jr-inc) > tmp_t )
    516                          t(jr) = t(jr-inc)
    517                          jr    = jr - inc
    518                          IF ( jr <= inc )  EXIT
    519                       ENDDO
    520                       t(jr) = tmp_t
    521                    ENDDO
    522                 ENDDO
    523 
    524          case3: DO  t_index = 1, t_index_number
    525 
    526                    pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
    527                    pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
    528                    pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
    529 
    530                    i3 = ( pos_x + 0.5_wp * dx ) * ddx
    531                    j3 = ( pos_y + 0.5_wp * dy ) * ddy
    532                    k3 = pos_z / dz + offset_ocean_nzt_m1
    533 
    534                    i5 = pos_x * ddx
    535                    j5 = pos_y * ddy
    536                    k5 = pos_z / dz + offset_ocean_nzt_m1
    537 
    538                    IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
    539                         nzb_s_inner(j5,i3) /= 0 )  THEN
    540 
    541                       IF ( pos_z == nzb_s_inner(j5,i3) * dz  .AND. &
    542                            k3 == nzb_s_inner(j5,i3) )  THEN
    543                          reflect_z = .TRUE.
    544                          EXIT case3
    545                       ENDIF
    546 
    547                    ENDIF
    548 
    549                    IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
    550                         nzb_s_inner(j3,i3) /= 0 )  THEN
    551 
    552                       IF ( pos_z == nzb_s_inner(j3,i3) * dz  .AND. &
    553                            k3 == nzb_s_inner(j3,i3) )  THEN
    554                          reflect_z = .TRUE.
    555                          EXIT case3
    556                       ENDIF
    557 
    558                       IF ( pos_y == ( j3 * dy - 0.5_wp * dy )  .AND. &
    559                            pos_z < nzb_s_inner(j3,i3) * dz )  THEN
    560                          reflect_y = .TRUE.
    561                          EXIT case3
    562                       ENDIF
    563 
    564                    ENDIF
    565 
    566                    IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
    567                         nzb_s_inner(j3,i5) /= 0 )  THEN
    568 
    569                       IF ( pos_z == nzb_s_inner(j3,i5) * dz  .AND. &
    570                            k3 == nzb_s_inner(j3,i5) )  THEN
    571                          reflect_z = .TRUE.
    572                          EXIT case3
    573                       ENDIF
    574 
    575                       IF ( pos_x == ( i5 * dx + 0.5_wp * dx )  .AND. &
    576                            pos_z < nzb_s_inner(j3,i5) * dz )  THEN
    577                          reflect_x = .TRUE.
    578                          EXIT case3
    579                       ENDIF
    580 
    581                    ENDIF
    582 
    583                 ENDDO case3
    584 
    585 !
    586 !--          Case 4
    587              ELSEIF ( particles(n)%x < pos_x_old  .AND. &
    588                       particles(n)%y < pos_y_old )  THEN
    589 
    590                 t_index = 1
    591 
    592                 DO  i = i1, i2, -1
    593                    xline      = i * dx - 0.5_wp * dx
    594                    t(t_index) = ( pos_x_old - xline ) / &
    595                                 ( pos_x_old - particles(n)%x )
    596                    t_index    = t_index + 1
    597                 ENDDO
    598 
    599                 DO  j = j1, j2, -1
    600                    yline      = j * dy - 0.5_wp * dy
    601                    t(t_index) = ( pos_y_old - yline ) / &
    602                                 ( pos_y_old - particles(n)%y )
    603                    t_index    = t_index + 1
    604                 ENDDO
    605 
    606                 IF ( particles(n)%z < pos_z_old )  THEN
    607                    DO  k = k1, k2 , -1
    608                       zline      = k * dz
    609                       t(t_index) = ( pos_z_old - zline ) / &
    610                                    ( pos_z_old-particles(n)%z )
    611                       t_index    = t_index + 1
    612                    ENDDO
    613                 ENDIF
    614                 t_index_number = t_index-1
    615 
    616 !
    617 !--             Sorting t
    618                 inc = 1
    619                 jr  = 1
    620 
    621                 DO WHILE ( inc <= t_index_number )
    622                    inc = 3 * inc + 1
    623                 ENDDO
    624 
    625                 DO WHILE ( inc > 1 )
    626                    inc = inc / 3
    627                    DO  ir = inc+1, t_index_number
    628                       tmp_t = t(ir)
    629                       jr = ir
    630                       DO WHILE ( t(jr-inc) > tmp_t )
    631                          t(jr) = t(jr-inc)
    632                          jr    = jr - inc
    633                          IF ( jr <= inc )  EXIT
    634                       ENDDO
    635                       t(jr) = tmp_t
    636                    ENDDO
    637                 ENDDO
    638 
    639          case4: DO  t_index = 1, t_index_number
    640 
    641                    pos_x = pos_x_old + t(t_index) * ( prt_x - pos_x_old )
    642                    pos_y = pos_y_old + t(t_index) * ( prt_y - pos_y_old )
    643                    pos_z = pos_z_old + t(t_index) * ( prt_z - pos_z_old )
    644 
    645                    i3 = ( pos_x + 0.5_wp * dx ) * ddx   
    646                    j3 = ( pos_y + 0.5_wp * dy ) * ddy
    647                    k3 = pos_z / dz + offset_ocean_nzt_m1
    648 
    649                    i5 = pos_x * ddx
    650                    j5 = pos_y * ddy
    651                    k5 = pos_z / dz + offset_ocean_nzt_m1
    652 
    653                    IF ( k3 <= nzb_s_inner(j3,i3)  .AND. &
    654                         nzb_s_inner(j3,i3) /= 0 )  THEN
    655 
    656                       IF ( pos_z == nzb_s_inner(j3,i3) * dz  .AND. &
    657                            k3 == nzb_s_inner(j3,i3) )  THEN
    658                          reflect_z = .TRUE.
    659                          EXIT case4
    660                       ENDIF
    661 
    662                    ENDIF
    663 
    664                    IF ( k5 <= nzb_s_inner(j3,i5)  .AND. &
    665                         nzb_s_inner(j3,i5) /= 0 )  THEN
    666 
    667                       IF ( pos_z == nzb_s_inner(j3,i5) * dz  .AND. &
    668                            k3 == nzb_s_inner(j3,i5) )  THEN
    669                          reflect_z = .TRUE.
    670                          EXIT case4
    671                       ENDIF
    672 
    673                       IF ( pos_x == ( i5 * dx + 0.5_wp * dx )  .AND. &
    674                            nzb_s_inner(j3,i5) /=0  .AND.          &
    675                            pos_z < nzb_s_inner(j3,i5) * dz )  THEN
    676                          reflect_x = .TRUE.
    677                          EXIT case4
    678                       ENDIF
    679 
    680                    ENDIF
    681 
    682                    IF ( k5 <= nzb_s_inner(j5,i3)  .AND. &
    683                         nzb_s_inner(j5,i3) /= 0 )  THEN
    684 
    685                       IF ( pos_z == nzb_s_inner(j5,i3) * dz  .AND. &
    686                            k5 == nzb_s_inner(j5,i3) )  THEN
    687                          reflect_z = .TRUE.
    688                          EXIT case4
    689                       ENDIF
    690 
    691                       IF ( pos_y == ( j5 * dy + 0.5_wp * dy )  .AND. &
    692                            nzb_s_inner(j5,i3) /= 0  .AND.         &
    693                            pos_z < nzb_s_inner(j5,i3) * dz )  THEN
    694                          reflect_y = .TRUE.
    695                          EXIT case4
    696                       ENDIF
    697 
    698                    ENDIF
    699  
    700                 ENDDO case4
     536!
     537!--             Swap time
     538                t_old = t(t_index)
     539
     540             ENDDO
     541!
     542!--          If a particle was reflected, calculate final position from last
     543!--          intermediate position.
     544             IF ( reflect_x .OR. reflect_y .OR. reflect_z )  THEN
     545
     546                particles(n)%x = pos_x + ( 1.0_wp - t_old ) * dt_particle      &
     547                                                         * particles(n)%speed_x
     548                particles(n)%y = pos_y + ( 1.0_wp - t_old ) * dt_particle      &
     549                                                         * particles(n)%speed_y
     550                particles(n)%z = pos_z + ( 1.0_wp - t_old ) * dt_particle      &
     551                                                         * particles(n)%speed_z
    701552
    702553             ENDIF
    703554
    704           ENDIF      ! Check, if particle position is below the surface
    705 
    706 !
    707 !--       Do the respective reflection, in case that one of the above conditions
    708 !--       is found to be true
    709           IF ( reflect_z )  THEN
    710 
    711              particles(n)%z       = 2.0_wp * pos_z - prt_z
    712              particles(n)%speed_z = - particles(n)%speed_z
    713 
    714              IF ( use_sgs_for_particles )  THEN
    715                 particles(n)%rvar3 = - particles(n)%rvar3
    716              ENDIF
    717              reflect_z = .FALSE.
    718 
    719           ELSEIF ( reflect_y )  THEN
    720 
    721              particles(n)%y       = 2.0_wp * pos_y - prt_y
    722              particles(n)%speed_y = - particles(n)%speed_y
    723 
    724              IF ( use_sgs_for_particles )  THEN
    725                 particles(n)%rvar2 = - particles(n)%rvar2
    726              ENDIF
    727              reflect_y = .FALSE.
    728 
    729           ELSEIF ( reflect_x )  THEN
    730 
    731              particles(n)%x       = 2.0_wp * pos_x - prt_x
    732              particles(n)%speed_x = - particles(n)%speed_x
    733 
    734              IF ( use_sgs_for_particles )  THEN
    735                 particles(n)%rvar1 = - particles(n)%rvar1
    736              ENDIF
    737              reflect_x = .FALSE.
    738 
    739           ENDIF       
    740    
     555          ENDIF
     556
    741557       ENDDO
    742558
Note: See TracChangeset for help on using the changeset viewer.