Changeset 2698 for palm/trunk


Ignore:
Timestamp:
Dec 14, 2017 6:46:24 PM (7 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

Location:
palm/trunk/SOURCE
Files:
14 edited

Legend:

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

    r2696 r2698  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Bugfix in definition of d_seconds_year.
    2323!
    2424! Former revisions:
     
    5656    REAL(wp) ::  time_utc_init = 43200.0_wp   !< UTC time at model start
    5757
    58     REAL(wp), PARAMETER ::  d_hours_day    = 0.0416666666667_wp    !< inverse of hours per day (1/24)
    59     REAL(wp), PARAMETER ::  d_seconds_hour = 0.000277777777778_wp   !< inverse of seconds per hour (1/3600)
    60     REAL(wp), PARAMETER ::  d_seconds_year = 31536000.0_wp          !< inverse of the seconds per year (1/(365*86400))
     58    REAL(wp), PARAMETER ::  d_hours_day    = 1.0_wp / 24.0_wp       !< inverse of hours per day (1/24)
     59    REAL(wp), PARAMETER ::  d_seconds_hour = 1.0_wp / 3600.0_wp     !< inverse of seconds per hour (1/3600)
     60    REAL(wp), PARAMETER ::  d_seconds_year = 1.0_wp / 31536000.0_wp !< inverse of the seconds per year (1/(365*86400))
    6161   
    6262    SAVE
  • palm/trunk/SOURCE/header.f90

    r2696 r2698  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Bugfix in get_topography_top_index
    2323!
    2424! Former revisions:
     
    410410
    411411    USE surface_mod,                                                           &
    412         ONLY:  surf_def_h, get_topography_top_index
     412        ONLY:  surf_def_h, get_topography_top_index_ji
    413413
    414414    USE synthetic_turbulence_generator_mod,                                    &
     
    993993       WRITE( io, 280 )
    994994       IF ( turbulent_inflow )  THEN
    995           WRITE( io, 281 )  zu( get_topography_top_index( 0, 0, 's' ) )
     995          WRITE( io, 281 )  zu( get_topography_top_index_ji( 0, 0, 's' ) )
    996996       ENDIF
    997997       IF ( TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
  • palm/trunk/SOURCE/init_3d_model.f90

    r2696 r2698  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Bugfix in get_topography_top_index
    2323!
    2424! Former revisions:
     
    466466    USE surface_mod,                                                           &
    467467        ONLY :  init_surface_arrays, init_surfaces, surf_def_h, surf_lsm_h,    &
    468                 surf_usm_h, get_topography_top_index
     468                surf_usm_h, get_topography_top_index_ji
    469469   
    470470    USE transpose_indices
     
    15801580          DO  i = nxlg, nxrg
    15811581             DO  j = nysg, nyng
    1582                 nz_u_shift = get_topography_top_index( j, i, 'u' )
    1583                 nz_v_shift = get_topography_top_index( j, i, 'v' )
    1584                 nz_w_shift = get_topography_top_index( j, i, 'w' )
    1585                 nz_s_shift = get_topography_top_index( j, i, 's' )
     1582                nz_u_shift = get_topography_top_index_ji( j, i, 'u' )
     1583                nz_v_shift = get_topography_top_index_ji( j, i, 'v' )
     1584                nz_w_shift = get_topography_top_index_ji( j, i, 'w' )
     1585                nz_s_shift = get_topography_top_index_ji( j, i, 's' )
    15861586
    15871587                u(nz_u_shift:nzt+1,j,i)  = u(0:nzt+1-nz_u_shift,j,i)               
     
    16251625          IF ( complex_terrain )  THEN
    16261626             IF ( nxlg <= 0 .AND. nxrg >= 0 .AND. nysg <= 0 .AND. nyng >= 0 )  THEN
    1627                 nz_u_shift_l = get_topography_top_index( 0, 0, 'u' )
    1628                 nz_v_shift_l = get_topography_top_index( 0, 0, 'v' )
    1629                 nz_w_shift_l = get_topography_top_index( 0, 0, 'w' )
    1630                 nz_s_shift_l = get_topography_top_index( 0, 0, 's' )
     1627                nz_u_shift_l = get_topography_top_index_ji( 0, 0, 'u' )
     1628                nz_v_shift_l = get_topography_top_index_ji( 0, 0, 'v' )
     1629                nz_w_shift_l = get_topography_top_index_ji( 0, 0, 'w' )
     1630                nz_s_shift_l = get_topography_top_index_ji( 0, 0, 's' )
    16311631             ELSE
    16321632                nz_u_shift_l = 0
  • palm/trunk/SOURCE/init_grid.f90

    r2696 r2698  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Bugfix in get_topography_top_index
    2323!
    2424! Former revisions:
     
    298298
    299299    USE surface_mod,                                                           &
    300         ONLY:  get_topography_top_index, init_bc
     300        ONLY:  get_topography_top_index, get_topography_top_index_ji, init_bc
    301301
    302302    USE vertical_nesting_mod,                                                  &
     
    552552!--          Topography height on scalar grid. Therefore, determine index of
    553553!--          upward-facing surface element on scalar grid.
    554              zu_s_inner(i,j) = zu( get_topography_top_index( j, i, 's' ) )
     554             zu_s_inner(i,j) = zu( get_topography_top_index_ji( j, i, 's' ) )
    555555!
    556556!--          Topography height on w grid. Therefore, determine index of
    557557!--          upward-facing surface element on w grid.
    558              zw_w_inner(i,j) = zw( get_topography_top_index( j, i, 's' ) )
     558             zw_w_inner(i,j) = zw( get_topography_top_index_ji( j, i, 's' ) )
    559559          ENDDO
    560560       ENDDO
     
    14711471
    14721472    USE surface_mod,                                                           &
    1473         ONLY:  get_topography_top_index
     1473        ONLY:  get_topography_top_index, get_topography_top_index_ji
    14741474
    14751475    IMPLICIT NONE
  • palm/trunk/SOURCE/lpm.f90

    r2696 r2698  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Grid indices passed to lpm_boundary_conds. (responsible Philipp Thiele)
    2323!
    2424! Former revisions:
     
    343343!--             optimization is still possible.)
    344344                IF ( topography /= 'flat' .AND. k < nzb_max + 2 )  THEN
    345                    CALL lpm_boundary_conds( 'walls' )
     345                   CALL  lpm_boundary_conds( 'walls', i, j, k )
    346346                ENDIF
    347347!
     
    353353!--             the top or bottom boundary and delete those particles, which are
    354354!--             older than allowed
    355                 CALL lpm_boundary_conds( 'bottom/top' )
     355                CALL lpm_boundary_conds( 'bottom/top', i, j, k )
    356356!
    357357!---            If not all particles of the actual grid cell have reached the
  • palm/trunk/SOURCE/lpm_advec.f90

    r2696 r2698  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Particle interpolations at walls in case of SGS velocities revised and not
     23! required parts are removed. (responsible Philipp Thiele)
     24!
     25! Bugfix in get_topography_top_index
    2326!
    2427! Former revisions:
     
    149152       
    150153    USE indices,                                                               &
    151         ONLY:  nzb, nzt
     154        ONLY:  nzb, nzt, wall_flags_0
    152155       
    153156    USE kinds
     
    163166
    164167    USE surface_mod,                                                           &
    165         ONLY:  get_topography_top_index, surf_def_h, surf_lsm_h, surf_usm_h
     168        ONLY:  get_topography_top_index_ji, surf_def_h, surf_lsm_h, surf_usm_h
    166169
    167170    IMPLICIT NONE
     171
     172    LOGICAL ::  subbox_at_wall !< flag to see if the current subgridbox is adjacent to a wall
    168173
    169174    INTEGER(iwp) ::  agp                         !< loop variable
     
    302307!--       above topography (Prandtl-layer height)
    303308!--       Determine vertical index of topography top
    304           k_wall = get_topography_top_index( jp, ip, 's' )
     309          k_wall = get_topography_top_index_ji( jp, ip, 's' )
    305310
    306311          IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
     
    397402!
    398403!--       Determine vertical index of topography top
    399           k_wall = get_topography_top_index( jp,ip, 's' )
     404          k_wall = get_topography_top_index_ji( jp,ip, 's' )
    400405
    401406          IF ( constant_flux_layer  .AND.  zv(n) - zw(k_wall) < z_p )  THEN
     
    528533!-- velocities
    529534    IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets )  THEN
    530 
    531        IF ( topography == 'flat' )  THEN
    532 
    533           DO  nb = 0,7
    534 
     535   
     536       DO  nb = 0,7
     537         
     538          subbox_at_wall = .FALSE.         
     539!
     540!--       In case of topography check if subbox is adjacent to a wall
     541          IF ( .NOT. topography == 'flat' ) THEN
     542             i = ip + MERGE( -1_iwp , 1_iwp, BTEST( nb, 2 ) )
     543             j = jp + MERGE( -1_iwp , 1_iwp, BTEST( nb, 1 ) )
     544             k = kp + MERGE( -1_iwp , 1_iwp, BTEST( nb, 0 ) )
     545             IF ( .NOT. BTEST(wall_flags_0(k,  jp, ip), 0) .OR.                &
     546                  .NOT. BTEST(wall_flags_0(kp, j,  ip), 0) .OR.                &
     547                  .NOT. BTEST(wall_flags_0(kp, jp, i ), 0) )                   &
     548             THEN
     549                subbox_at_wall = .TRUE.
     550             ENDIF
     551          ENDIF
     552          IF ( subbox_at_wall ) THEN
     553             e_int(start_index(nb):end_index(nb))     = e(kp,jp,ip)
     554             diss_int(start_index(nb):end_index(nb))  = diss(kp,jp,ip)
     555             de_dx_int(start_index(nb):end_index(nb)) = de_dx(kp,jp,ip)
     556             de_dy_int(start_index(nb):end_index(nb)) = de_dy(kp,jp,ip)
     557             de_dz_int(start_index(nb):end_index(nb)) = de_dz(kp,jp,ip)
     558!
     559!--          Set flag for stochastic equation.
     560             term_1_2(start_index(nb):end_index(nb)) = 0.0_wp             
     561          ELSE
    535562             i = ip + block_offset(nb)%i_off
    536563             j = jp + block_offset(nb)%j_off
     
    600627                ELSE
    601628                   de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
    602                                   ( gg - bb ) * de_dy(k+1,j,i+1) + &
    603                                   ( gg - cc ) * de_dy(k+1,j+1,i) + &
    604                                   ( gg - dd ) * de_dy(k+1,j+1,i+1) &
     629                                   ( gg - bb ) * de_dy(k+1,j,i+1) + &
     630                                   ( gg - cc ) * de_dy(k+1,j+1,i) + &
     631                                   ( gg - dd ) * de_dy(k+1,j+1,i+1) &
    605632                                  ) / ( 3.0_wp * gg )
    606                    de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / dz * &
    607                                               ( de_dy_int_u - de_dy_int_l )
     633                      de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / dz * &
     634                                                 ( de_dy_int_u - de_dy_int_l )
    608635                ENDIF
    609636
     
    638665                               ( gg - cc ) * diss(k,j+1,i) + &
    639666                               ( gg - dd ) * diss(k,j+1,i+1) &
    640                               ) / ( 3.0_wp * gg )
     667                               ) / ( 3.0_wp * gg )
    641668
    642669                IF ( k == nzt )  THEN
     
    649676                                 ) / ( 3.0_wp * gg )
    650677                   diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz * &
    651                                            ( diss_int_u - diss_int_l )
     678                                            ( diss_int_u - diss_int_l )
    652679                ENDIF
    653680
     
    655682!--             Set flag for stochastic equation.
    656683                term_1_2(n) = 1.0_wp
    657 
    658684             ENDDO
    659           ENDDO
    660 
    661        ELSE ! non-flat topography, e.g., buildings
    662 
    663           DO  nb = 0, 7
    664 
    665              DO  n = start_index(nb), end_index(nb)
    666 
    667                 i = particles(n)%x * ddx
    668                 j = particles(n)%y * ddy
    669                 k = ( zv(n) + dz * atmos_ocean_sign ) / dz  &
    670                     + offset_ocean_nzt                      ! only exact if eq.dist
    671 !
    672 !--             In case that there are buildings it has to be determined
    673 !--             how many of the gridpoints defining the particle box are
    674 !--             situated within a building
    675 !--             gp_outside_of_building(1): i,j,k
    676 !--             gp_outside_of_building(2): i,j+1,k
    677 !--             gp_outside_of_building(3): i,j,k+1
    678 !--             gp_outside_of_building(4): i,j+1,k+1
    679 !--             gp_outside_of_building(5): i+1,j,k
    680 !--             gp_outside_of_building(6): i+1,j+1,k
    681 !--             gp_outside_of_building(7): i+1,j,k+1
    682 !--             gp_outside_of_building(8): i+1,j+1,k+1
    683 
    684                 gp_outside_of_building = 0
    685                 location = 0.0_wp
    686                 num_gp = 0
    687 
    688 !
    689 !--             Determine vertical index of topography top at (j,i)
    690                 k_wall = get_topography_top_index( j, i, 's' )
    691 !
    692 !--             To do: Reconsider order of computations in order to avoid
    693 !--             unnecessary index calculations.
    694                 IF ( k > k_wall  .OR.  k_wall == 0 )  THEN
    695                    num_gp = num_gp + 1
    696                    gp_outside_of_building(1) = 1
    697                    location(num_gp,1) = i * dx
    698                    location(num_gp,2) = j * dy
    699                    location(num_gp,3) = k * dz - 0.5_wp * dz
    700                    ei(num_gp)     = e(k,j,i)
    701                    dissi(num_gp)  = diss(k,j,i)
    702                    de_dxi(num_gp) = de_dx(k,j,i)
    703                    de_dyi(num_gp) = de_dy(k,j,i)
    704                    de_dzi(num_gp) = de_dz(k,j,i)
    705                 ENDIF
    706 
    707 !
    708 !--             Determine vertical index of topography top at (j+1,i)
    709                 k_wall = get_topography_top_index( j+1, i, 's' )
    710 
    711                 IF ( k > k_wall  .OR.  k_wall == 0 )  THEN
    712                    num_gp = num_gp + 1
    713                    gp_outside_of_building(2) = 1
    714                    location(num_gp,1) = i * dx
    715                    location(num_gp,2) = (j+1) * dy
    716                    location(num_gp,3) = k * dz - 0.5_wp * dz
    717                    ei(num_gp)     = e(k,j+1,i)
    718                    dissi(num_gp)  = diss(k,j+1,i)
    719                    de_dxi(num_gp) = de_dx(k,j+1,i)
    720                    de_dyi(num_gp) = de_dy(k,j+1,i)
    721                    de_dzi(num_gp) = de_dz(k,j+1,i)
    722                 ENDIF
    723 
    724 !
    725 !--             Determine vertical index of topography top at (j,i)
    726                 k_wall = get_topography_top_index( j, i, 's' )
    727 
    728                 IF ( k+1 > k_wall  .OR.  k_wall == 0 )  THEN
    729                    num_gp = num_gp + 1
    730                    gp_outside_of_building(3) = 1
    731                    location(num_gp,1) = i * dx
    732                    location(num_gp,2) = j * dy
    733                    location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
    734                    ei(num_gp)     = e(k+1,j,i)
    735                    dissi(num_gp)  = diss(k+1,j,i)
    736                    de_dxi(num_gp) = de_dx(k+1,j,i)
    737                    de_dyi(num_gp) = de_dy(k+1,j,i)
    738                    de_dzi(num_gp) = de_dz(k+1,j,i)
    739                 ENDIF
    740 
    741 !
    742 !--             Determine vertical index of topography top at (j+1,i)
    743                 k_wall = get_topography_top_index( j+1, i, 's' )
    744                 IF ( k+1 > k_wall  .OR.  k_wall == 0 )  THEN
    745                    num_gp = num_gp + 1
    746                    gp_outside_of_building(4) = 1
    747                    location(num_gp,1) = i * dx
    748                    location(num_gp,2) = (j+1) * dy
    749                    location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
    750                    ei(num_gp)     = e(k+1,j+1,i)
    751                    dissi(num_gp)  = diss(k+1,j+1,i)
    752                    de_dxi(num_gp) = de_dx(k+1,j+1,i)
    753                    de_dyi(num_gp) = de_dy(k+1,j+1,i)
    754                    de_dzi(num_gp) = de_dz(k+1,j+1,i)
    755                 ENDIF
    756 
    757 !
    758 !--             Determine vertical index of topography top at (j,i+1)
    759                 k_wall = get_topography_top_index( j, i+1, 's' )
    760                 IF ( k > k_wall  .OR.  k_wall == 0 )  THEN
    761                    num_gp = num_gp + 1
    762                    gp_outside_of_building(5) = 1
    763                    location(num_gp,1) = (i+1) * dx
    764                    location(num_gp,2) = j * dy
    765                    location(num_gp,3) = k * dz - 0.5_wp * dz
    766                    ei(num_gp)     = e(k,j,i+1)
    767                    dissi(num_gp)  = diss(k,j,i+1)
    768                    de_dxi(num_gp) = de_dx(k,j,i+1)
    769                    de_dyi(num_gp) = de_dy(k,j,i+1)
    770                    de_dzi(num_gp) = de_dz(k,j,i+1)
    771                 ENDIF
    772 
    773 !
    774 !--             Determine vertical index of topography top at (j+1,i+1)
    775                 k_wall = get_topography_top_index( j+1, i+1, 's' )
    776 
    777                 IF ( k > k_wall  .OR.  k_wall == 0 )  THEN
    778                    num_gp = num_gp + 1
    779                    gp_outside_of_building(6) = 1
    780                    location(num_gp,1) = (i+1) * dx
    781                    location(num_gp,2) = (j+1) * dy
    782                    location(num_gp,3) = k * dz - 0.5_wp * dz
    783                    ei(num_gp)     = e(k,j+1,i+1)
    784                    dissi(num_gp)  = diss(k,j+1,i+1)
    785                    de_dxi(num_gp) = de_dx(k,j+1,i+1)
    786                    de_dyi(num_gp) = de_dy(k,j+1,i+1)
    787                    de_dzi(num_gp) = de_dz(k,j+1,i+1)
    788                 ENDIF
    789 
    790 !
    791 !--             Determine vertical index of topography top at (j,i+1)
    792                 k_wall = get_topography_top_index( j, i+1, 's' )
    793 
    794                 IF ( k+1 > k_wall  .OR.  k_wall == 0 )  THEN
    795                    num_gp = num_gp + 1
    796                    gp_outside_of_building(7) = 1
    797                    location(num_gp,1) = (i+1) * dx
    798                    location(num_gp,2) = j * dy
    799                    location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
    800                    ei(num_gp)     = e(k+1,j,i+1)
    801                    dissi(num_gp)  = diss(k+1,j,i+1)
    802                    de_dxi(num_gp) = de_dx(k+1,j,i+1)
    803                    de_dyi(num_gp) = de_dy(k+1,j,i+1)
    804                    de_dzi(num_gp) = de_dz(k+1,j,i+1)
    805                 ENDIF
    806 
    807 !
    808 !--             Determine vertical index of topography top at (j+1,i+1)
    809                 k_wall = get_topography_top_index( j+1, i+1, 's' )
    810 
    811                 IF ( k+1 > k_wall  .OR.  k_wall == 0)  THEN
    812                    num_gp = num_gp + 1
    813                    gp_outside_of_building(8) = 1
    814                    location(num_gp,1) = (i+1) * dx
    815                    location(num_gp,2) = (j+1) * dy
    816                    location(num_gp,3) = (k+1) * dz - 0.5_wp * dz
    817                    ei(num_gp)     = e(k+1,j+1,i+1)
    818                    dissi(num_gp)  = diss(k+1,j+1,i+1)
    819                    de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
    820                    de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
    821                    de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
    822                 ENDIF
    823 !
    824 !--             If all gridpoints are situated outside of a building, then the
    825 !--             ordinary interpolation scheme can be used.
    826                 IF ( num_gp == 8 )  THEN
    827 
    828                    x  = particles(n)%x - i * dx
    829                    y  = particles(n)%y - j * dy
    830                    aa = x**2          + y**2
    831                    bb = ( dx - x )**2 + y**2
    832                    cc = x**2          + ( dy - y )**2
    833                    dd = ( dx - x )**2 + ( dy - y )**2
    834                    gg = aa + bb + cc + dd
    835      
    836                    e_int_l = ( ( gg - aa ) * e(k,j,i)   + ( gg - bb ) * e(k,j,i+1)   &
    837                              + ( gg - cc ) * e(k,j+1,i) + ( gg - dd ) * e(k,j+1,i+1) &
    838                              ) / ( 3.0_wp * gg )
    839      
    840                    IF ( k == nzt )  THEN
    841                       e_int(n) = e_int_l
    842                    ELSE
    843                       e_int_u = ( ( gg - aa ) * e(k+1,j,i)   + &
    844                                   ( gg - bb ) * e(k+1,j,i+1) + &
    845                                   ( gg - cc ) * e(k+1,j+1,i) + &
    846                                   ( gg - dd ) * e(k+1,j+1,i+1) &
    847                                 ) / ( 3.0_wp * gg )
    848                       e_int(n) = e_int_l + ( zv(n) - zu(k) ) / dz * &
    849                                           ( e_int_u - e_int_l )
    850                    ENDIF
    851 !
    852 !--                Needed to avoid NaN particle velocities (this might not be
    853 !--                required any more)
    854                    IF ( e_int(n) <= 0.0_wp )  THEN
    855                       e_int(n) = 1.0E-20_wp
    856                    ENDIF
    857 !
    858 !--                Interpolate the TKE gradient along x (adopt incides i,j,k
    859 !--                and all position variables from above (TKE))
    860                    de_dx_int_l = ( ( gg - aa ) * de_dx(k,j,i)   + &
    861                                    ( gg - bb ) * de_dx(k,j,i+1) + &
    862                                    ( gg - cc ) * de_dx(k,j+1,i) + &
    863                                    ( gg - dd ) * de_dx(k,j+1,i+1) &
    864                                  ) / ( 3.0_wp * gg )
    865 
    866                    IF ( ( k == nzt )  .OR.  ( k == nzb ) )  THEN
    867                       de_dx_int(n) = de_dx_int_l
    868                    ELSE
    869                       de_dx_int_u = ( ( gg - aa ) * de_dx(k+1,j,i)   + &
    870                                       ( gg - bb ) * de_dx(k+1,j,i+1) + &
    871                                       ( gg - cc ) * de_dx(k+1,j+1,i) + &
    872                                       ( gg - dd ) * de_dx(k+1,j+1,i+1) &
    873                                     ) / ( 3.0_wp * gg )
    874                       de_dx_int(n) = de_dx_int_l + ( zv(n) - zu(k) ) / &
    875                                               dz * ( de_dx_int_u - de_dx_int_l )
    876                    ENDIF
    877 
    878 !
    879 !--                Interpolate the TKE gradient along y
    880                    de_dy_int_l = ( ( gg - aa ) * de_dy(k,j,i)   + &
    881                                    ( gg - bb ) * de_dy(k,j,i+1) + &
    882                                    ( gg - cc ) * de_dy(k,j+1,i) + &
    883                                    ( gg - dd ) * de_dy(k,j+1,i+1) &
    884                                  ) / ( 3.0_wp * gg )
    885                    IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
    886                       de_dy_int(n) = de_dy_int_l
    887                    ELSE
    888                       de_dy_int_u = ( ( gg - aa ) * de_dy(k+1,j,i)   + &
    889                                       ( gg - bb ) * de_dy(k+1,j,i+1) + &
    890                                       ( gg - cc ) * de_dy(k+1,j+1,i) + &
    891                                       ( gg - dd ) * de_dy(k+1,j+1,i+1) &
    892                                     ) / ( 3.0_wp * gg )
    893                       de_dy_int(n) = de_dy_int_l + ( zv(n) - zu(k) ) / &
    894                                               dz * ( de_dy_int_u - de_dy_int_l )
    895                    ENDIF
    896 
    897 !
    898 !--                Interpolate the TKE gradient along z
    899                    IF ( zv(n) < 0.5_wp * dz )  THEN
    900                       de_dz_int(n) = 0.0_wp
    901                    ELSE
    902                       de_dz_int_l = ( ( gg - aa ) * de_dz(k,j,i)   + &
    903                                       ( gg - bb ) * de_dz(k,j,i+1) + &
    904                                       ( gg - cc ) * de_dz(k,j+1,i) + &
    905                                       ( gg - dd ) * de_dz(k,j+1,i+1) &
    906                                     ) / ( 3.0_wp * gg )
    907    
    908                       IF ( ( k+1 == nzt+1 )  .OR.  ( k == nzb ) )  THEN
    909                          de_dz_int(n) = de_dz_int_l
    910                       ELSE
    911                          de_dz_int_u = ( ( gg - aa ) * de_dz(k+1,j,i)   + &
    912                                          ( gg - bb ) * de_dz(k+1,j,i+1) + &
    913                                          ( gg - cc ) * de_dz(k+1,j+1,i) + &
    914                                          ( gg - dd ) * de_dz(k+1,j+1,i+1) &
    915                                        ) / ( 3.0_wp * gg )
    916                          de_dz_int(n) = de_dz_int_l + ( zv(n) - zu(k) ) /&
    917                                               dz * ( de_dz_int_u - de_dz_int_l )
    918                       ENDIF
    919                    ENDIF
    920 
    921 !
    922 !--                Interpolate the dissipation of TKE
    923                    diss_int_l = ( ( gg - aa ) * diss(k,j,i)   + &
    924                                   ( gg - bb ) * diss(k,j,i+1) + &
    925                                   ( gg - cc ) * diss(k,j+1,i) + &
    926                                   ( gg - dd ) * diss(k,j+1,i+1) &
    927                                 ) / ( 3.0_wp * gg )
    928    
    929                    IF ( k == nzt )  THEN
    930                       diss_int(n) = diss_int_l
    931                    ELSE
    932                       diss_int_u = ( ( gg - aa ) * diss(k+1,j,i)   + &
    933                                      ( gg - bb ) * diss(k+1,j,i+1) + &
    934                                      ( gg - cc ) * diss(k+1,j+1,i) + &
    935                                      ( gg - dd ) * diss(k+1,j+1,i+1) &
    936                                    ) / ( 3.0_wp * gg )
    937                       diss_int(n) = diss_int_l + ( zv(n) - zu(k) ) / dz *&
    938                                               ( diss_int_u - diss_int_l )
    939                    ENDIF
    940 !
    941 !--                Set flag for stochastic equation.
    942                    term_1_2(n) = 1.0_wp
    943    
    944                 ELSE
    945    
    946 !
    947 !--                If wall between gridpoint 1 and gridpoint 5, then
    948 !--                Neumann boundary condition has to be applied
    949                    IF ( gp_outside_of_building(1) == 1  .AND. &
    950                         gp_outside_of_building(5) == 0 )  THEN
    951                       num_gp = num_gp + 1
    952                       location(num_gp,1) = i * dx + 0.5_wp * dx
    953                       location(num_gp,2) = j * dy
    954                       location(num_gp,3) = k * dz - 0.5_wp * dz
    955                       ei(num_gp)     = e(k,j,i)
    956                       dissi(num_gp)  = diss(k,j,i)
    957                       de_dxi(num_gp) = 0.0_wp
    958                       de_dyi(num_gp) = de_dy(k,j,i)
    959                       de_dzi(num_gp) = de_dz(k,j,i)
    960                    ENDIF
    961 
    962                    IF ( gp_outside_of_building(5) == 1  .AND. &
    963                         gp_outside_of_building(1) == 0 )  THEN
    964                       num_gp = num_gp + 1
    965                       location(num_gp,1) = i * dx + 0.5_wp * dx
    966                       location(num_gp,2) = j * dy
    967                       location(num_gp,3) = k * dz - 0.5_wp * dz
    968                       ei(num_gp)     = e(k,j,i+1)
    969                       dissi(num_gp)  = diss(k,j,i+1)
    970                       de_dxi(num_gp) = 0.0_wp
    971                       de_dyi(num_gp) = de_dy(k,j,i+1)
    972                       de_dzi(num_gp) = de_dz(k,j,i+1)
    973                    ENDIF
    974 
    975 !
    976 !--                If wall between gridpoint 5 and gridpoint 6, then
    977 !--                then Neumann boundary condition has to be applied
    978                    IF ( gp_outside_of_building(5) == 1  .AND. &
    979                         gp_outside_of_building(6) == 0 )  THEN
    980                       num_gp = num_gp + 1
    981                       location(num_gp,1) = (i+1) * dx
    982                       location(num_gp,2) = j * dy + 0.5_wp * dy
    983                       location(num_gp,3) = k * dz - 0.5_wp * dz
    984                       ei(num_gp)     = e(k,j,i+1)
    985                       dissi(num_gp)  = diss(k,j,i+1)
    986                       de_dxi(num_gp) = de_dx(k,j,i+1)
    987                       de_dyi(num_gp) = 0.0_wp
    988                       de_dzi(num_gp) = de_dz(k,j,i+1)
    989                    ENDIF
    990 
    991                    IF ( gp_outside_of_building(6) == 1  .AND. &
    992                         gp_outside_of_building(5) == 0 )  THEN
    993                       num_gp = num_gp + 1
    994                       location(num_gp,1) = (i+1) * dx
    995                       location(num_gp,2) = j * dy + 0.5_wp * dy
    996                       location(num_gp,3) = k * dz - 0.5_wp * dz
    997                       ei(num_gp)     = e(k,j+1,i+1)
    998                       dissi(num_gp)  = diss(k,j+1,i+1)
    999                       de_dxi(num_gp) = de_dx(k,j+1,i+1)
    1000                       de_dyi(num_gp) = 0.0_wp
    1001                       de_dzi(num_gp) = de_dz(k,j+1,i+1)
    1002                    ENDIF
    1003 
    1004 !
    1005 !--                If wall between gridpoint 2 and gridpoint 6, then
    1006 !--                Neumann boundary condition has to be applied
    1007                    IF ( gp_outside_of_building(2) == 1  .AND. &
    1008                         gp_outside_of_building(6) == 0 )  THEN
    1009                       num_gp = num_gp + 1
    1010                       location(num_gp,1) = i * dx + 0.5_wp * dx
    1011                       location(num_gp,2) = (j+1) * dy
    1012                       location(num_gp,3) = k * dz - 0.5_wp * dz
    1013                       ei(num_gp)     = e(k,j+1,i)
    1014                       dissi(num_gp)  = diss(k,j+1,i)
    1015                       de_dxi(num_gp) = 0.0_wp
    1016                       de_dyi(num_gp) = de_dy(k,j+1,i)
    1017                       de_dzi(num_gp) = de_dz(k,j+1,i)
    1018                    ENDIF
    1019 
    1020                    IF ( gp_outside_of_building(6) == 1  .AND. &
    1021                         gp_outside_of_building(2) == 0 )  THEN
    1022                       num_gp = num_gp + 1
    1023                       location(num_gp,1) = i * dx + 0.5_wp * dx
    1024                       location(num_gp,2) = (j+1) * dy
    1025                       location(num_gp,3) = k * dz - 0.5_wp * dz
    1026                       ei(num_gp)     = e(k,j+1,i+1)
    1027                       dissi(num_gp)  = diss(k,j+1,i+1)
    1028                       de_dxi(num_gp) = 0.0_wp
    1029                       de_dyi(num_gp) = de_dy(k,j+1,i+1)
    1030                       de_dzi(num_gp) = de_dz(k,j+1,i+1)
    1031                    ENDIF
    1032 
    1033 !
    1034 !--                If wall between gridpoint 1 and gridpoint 2, then
    1035 !--                Neumann boundary condition has to be applied
    1036                    IF ( gp_outside_of_building(1) == 1  .AND. &
    1037                         gp_outside_of_building(2) == 0 )  THEN
    1038                       num_gp = num_gp + 1
    1039                       location(num_gp,1) = i * dx
    1040                       location(num_gp,2) = j * dy + 0.5_wp * dy
    1041                       location(num_gp,3) = k * dz - 0.5_wp * dz
    1042                       ei(num_gp)     = e(k,j,i)
    1043                       dissi(num_gp)  = diss(k,j,i)
    1044                       de_dxi(num_gp) = de_dx(k,j,i)
    1045                       de_dyi(num_gp) = 0.0_wp
    1046                       de_dzi(num_gp) = de_dz(k,j,i) 
    1047                    ENDIF
    1048 
    1049                    IF ( gp_outside_of_building(2) == 1  .AND. &
    1050                         gp_outside_of_building(1) == 0 )  THEN
    1051                       num_gp = num_gp + 1
    1052                       location(num_gp,1) = i * dx
    1053                       location(num_gp,2) = j * dy + 0.5_wp * dy
    1054                       location(num_gp,3) = k * dz - 0.5_wp * dz
    1055                       ei(num_gp)     = e(k,j+1,i)
    1056                       dissi(num_gp)  = diss(k,j+1,i)
    1057                       de_dxi(num_gp) = de_dx(k,j+1,i)
    1058                       de_dyi(num_gp) = 0.0_wp
    1059                       de_dzi(num_gp) = de_dz(k,j+1,i)
    1060                    ENDIF
    1061 
    1062 !
    1063 !--                If wall between gridpoint 3 and gridpoint 7, then
    1064 !--                Neumann boundary condition has to be applied
    1065                    IF ( gp_outside_of_building(3) == 1  .AND. &
    1066                         gp_outside_of_building(7) == 0 )  THEN
    1067                       num_gp = num_gp + 1
    1068                       location(num_gp,1) = i * dx + 0.5_wp * dx
    1069                       location(num_gp,2) = j * dy
    1070                       location(num_gp,3) = k * dz + 0.5_wp * dz
    1071                       ei(num_gp)     = e(k+1,j,i)
    1072                       dissi(num_gp)  = diss(k+1,j,i)
    1073                       de_dxi(num_gp) = 0.0_wp
    1074                       de_dyi(num_gp) = de_dy(k+1,j,i)
    1075                       de_dzi(num_gp) = de_dz(k+1,j,i)
    1076                    ENDIF
    1077 
    1078                    IF ( gp_outside_of_building(7) == 1  .AND. &
    1079                         gp_outside_of_building(3) == 0 )  THEN
    1080                       num_gp = num_gp + 1
    1081                       location(num_gp,1) = i * dx + 0.5_wp * dx
    1082                       location(num_gp,2) = j * dy
    1083                       location(num_gp,3) = k * dz + 0.5_wp * dz
    1084                       ei(num_gp)     = e(k+1,j,i+1)
    1085                       dissi(num_gp)  = diss(k+1,j,i+1)
    1086                       de_dxi(num_gp) = 0.0_wp
    1087                       de_dyi(num_gp) = de_dy(k+1,j,i+1)
    1088                       de_dzi(num_gp) = de_dz(k+1,j,i+1)
    1089                    ENDIF
    1090 
    1091 !
    1092 !--                If wall between gridpoint 7 and gridpoint 8, then
    1093 !--                Neumann boundary condition has to be applied
    1094                    IF ( gp_outside_of_building(7) == 1  .AND. &
    1095                         gp_outside_of_building(8) == 0 )  THEN
    1096                       num_gp = num_gp + 1
    1097                       location(num_gp,1) = (i+1) * dx
    1098                       location(num_gp,2) = j * dy + 0.5_wp * dy
    1099                       location(num_gp,3) = k * dz + 0.5_wp * dz
    1100                       ei(num_gp)     = e(k+1,j,i+1)
    1101                       dissi(num_gp)  = diss(k+1,j,i+1)
    1102                       de_dxi(num_gp) = de_dx(k+1,j,i+1)
    1103                       de_dyi(num_gp) = 0.0_wp
    1104                       de_dzi(num_gp) = de_dz(k+1,j,i+1)
    1105                    ENDIF
    1106 
    1107                    IF ( gp_outside_of_building(8) == 1  .AND. &
    1108                         gp_outside_of_building(7) == 0 )  THEN
    1109                       num_gp = num_gp + 1
    1110                       location(num_gp,1) = (i+1) * dx
    1111                       location(num_gp,2) = j * dy + 0.5_wp * dy
    1112                       location(num_gp,3) = k * dz + 0.5_wp * dz
    1113                       ei(num_gp)     = e(k+1,j+1,i+1)
    1114                       dissi(num_gp)  = diss(k+1,j+1,i+1)
    1115                       de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
    1116                       de_dyi(num_gp) = 0.0_wp
    1117                       de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
    1118                    ENDIF
    1119 
    1120 !
    1121 !--                If wall between gridpoint 4 and gridpoint 8, then
    1122 !--                Neumann boundary condition has to be applied
    1123                    IF ( gp_outside_of_building(4) == 1  .AND. &
    1124                         gp_outside_of_building(8) == 0 )  THEN
    1125                       num_gp = num_gp + 1
    1126                       location(num_gp,1) = i * dx + 0.5_wp * dx
    1127                       location(num_gp,2) = (j+1) * dy
    1128                       location(num_gp,3) = k * dz + 0.5_wp * dz
    1129                       ei(num_gp)     = e(k+1,j+1,i)
    1130                       dissi(num_gp)  = diss(k+1,j+1,i)
    1131                       de_dxi(num_gp) = 0.0_wp
    1132                       de_dyi(num_gp) = de_dy(k+1,j+1,i)
    1133                       de_dzi(num_gp) = de_dz(k+1,j+1,i)
    1134                    ENDIF
    1135 
    1136                    IF ( gp_outside_of_building(8) == 1  .AND. &
    1137                         gp_outside_of_building(4) == 0 )  THEN
    1138                       num_gp = num_gp + 1
    1139                       location(num_gp,1) = i * dx + 0.5_wp * dx
    1140                       location(num_gp,2) = (j+1) * dy
    1141                       location(num_gp,3) = k * dz + 0.5_wp * dz
    1142                       ei(num_gp)     = e(k+1,j+1,i+1)
    1143                       dissi(num_gp)  = diss(k+1,j+1,i+1)
    1144                       de_dxi(num_gp) = 0.0_wp
    1145                       de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
    1146                       de_dzi(num_gp) = de_dz(k+1,j+1,i+1)
    1147                    ENDIF
    1148 
    1149 !
    1150 !--                If wall between gridpoint 3 and gridpoint 4, then
    1151 !--                Neumann boundary condition has to be applied
    1152                    IF ( gp_outside_of_building(3) == 1  .AND. &
    1153                         gp_outside_of_building(4) == 0 )  THEN
    1154                       num_gp = num_gp + 1
    1155                       location(num_gp,1) = i * dx
    1156                       location(num_gp,2) = j * dy + 0.5_wp * dy
    1157                       location(num_gp,3) = k * dz + 0.5_wp * dz
    1158                       ei(num_gp)     = e(k+1,j,i)
    1159                       dissi(num_gp)  = diss(k+1,j,i)
    1160                       de_dxi(num_gp) = de_dx(k+1,j,i)
    1161                       de_dyi(num_gp) = 0.0_wp
    1162                       de_dzi(num_gp) = de_dz(k+1,j,i)
    1163                    ENDIF
    1164 
    1165                    IF ( gp_outside_of_building(4) == 1  .AND. &
    1166                         gp_outside_of_building(3) == 0 )  THEN
    1167                       num_gp = num_gp + 1
    1168                       location(num_gp,1) = i * dx
    1169                       location(num_gp,2) = j * dy + 0.5_wp * dy
    1170                       location(num_gp,3) = k * dz + 0.5_wp * dz
    1171                       ei(num_gp)     = e(k+1,j+1,i)
    1172                       dissi(num_gp)  = diss(k+1,j+1,i)
    1173                       de_dxi(num_gp) = de_dx(k+1,j+1,i)
    1174                       de_dyi(num_gp) = 0.0_wp
    1175                       de_dzi(num_gp) = de_dz(k+1,j+1,i)
    1176                    ENDIF
    1177 
    1178 !
    1179 !--                If wall between gridpoint 1 and gridpoint 3, then
    1180 !--                Neumann boundary condition has to be applied
    1181 !--                (only one case as only building beneath is possible)
    1182                    IF ( gp_outside_of_building(1) == 0  .AND. &
    1183                         gp_outside_of_building(3) == 1 )  THEN
    1184                       num_gp = num_gp + 1
    1185                       location(num_gp,1) = i * dx
    1186                       location(num_gp,2) = j * dy
    1187                       location(num_gp,3) = k * dz
    1188                       ei(num_gp)     = e(k+1,j,i)
    1189                       dissi(num_gp)  = diss(k+1,j,i)
    1190                       de_dxi(num_gp) = de_dx(k+1,j,i)
    1191                       de_dyi(num_gp) = de_dy(k+1,j,i)
    1192                       de_dzi(num_gp) = 0.0_wp
    1193                    ENDIF
    1194 
    1195 !
    1196 !--                If wall between gridpoint 5 and gridpoint 7, then
    1197 !--                Neumann boundary condition has to be applied
    1198 !--                (only one case as only building beneath is possible)
    1199                    IF ( gp_outside_of_building(5) == 0  .AND. &
    1200                         gp_outside_of_building(7) == 1 )  THEN
    1201                       num_gp = num_gp + 1
    1202                       location(num_gp,1) = (i+1) * dx
    1203                       location(num_gp,2) = j * dy
    1204                       location(num_gp,3) = k * dz
    1205                       ei(num_gp)     = e(k+1,j,i+1)
    1206                       dissi(num_gp)  = diss(k+1,j,i+1)
    1207                       de_dxi(num_gp) = de_dx(k+1,j,i+1)
    1208                       de_dyi(num_gp) = de_dy(k+1,j,i+1)
    1209                       de_dzi(num_gp) = 0.0_wp
    1210                    ENDIF
    1211 
    1212 !
    1213 !--                If wall between gridpoint 2 and gridpoint 4, then
    1214 !--                Neumann boundary condition has to be applied
    1215 !--                (only one case as only building beneath is possible)
    1216                    IF ( gp_outside_of_building(2) == 0  .AND. &
    1217                         gp_outside_of_building(4) == 1 )  THEN
    1218                       num_gp = num_gp + 1
    1219                       location(num_gp,1) = i * dx
    1220                       location(num_gp,2) = (j+1) * dy
    1221                       location(num_gp,3) = k * dz
    1222                       ei(num_gp)     = e(k+1,j+1,i)
    1223                       dissi(num_gp)  = diss(k+1,j+1,i)
    1224                       de_dxi(num_gp) = de_dx(k+1,j+1,i)
    1225                       de_dyi(num_gp) = de_dy(k+1,j+1,i)
    1226                       de_dzi(num_gp) = 0.0_wp
    1227                    ENDIF
    1228 
    1229 !
    1230 !--                If wall between gridpoint 6 and gridpoint 8, then
    1231 !--                Neumann boundary condition has to be applied
    1232 !--                (only one case as only building beneath is possible)
    1233                    IF ( gp_outside_of_building(6) == 0  .AND. &
    1234                         gp_outside_of_building(8) == 1 )  THEN
    1235                       num_gp = num_gp + 1
    1236                       location(num_gp,1) = (i+1) * dx
    1237                       location(num_gp,2) = (j+1) * dy
    1238                       location(num_gp,3) = k * dz
    1239                       ei(num_gp)     = e(k+1,j+1,i+1)
    1240                       dissi(num_gp)  = diss(k+1,j+1,i+1)
    1241                       de_dxi(num_gp) = de_dx(k+1,j+1,i+1)
    1242                       de_dyi(num_gp) = de_dy(k+1,j+1,i+1)
    1243                       de_dzi(num_gp) = 0.0_wp
    1244                    ENDIF
    1245      
    1246 !
    1247 !--                Carry out the interpolation
    1248                    IF ( num_gp == 1 )  THEN
    1249 !
    1250 !--                   If only one of the gridpoints is situated outside of the
    1251 !--                   building, it follows that the values at the particle
    1252 !--                   location are the same as the gridpoint values
    1253                       e_int(n)     = ei(num_gp)   
    1254                       diss_int(n)  = dissi(num_gp)
    1255                       de_dx_int(n) = de_dxi(num_gp)
    1256                       de_dy_int(n) = de_dyi(num_gp)
    1257                       de_dz_int(n) = de_dzi(num_gp)
    1258 !
    1259 !--                   Set flag for stochastic equation which disables calculation
    1260 !--                   of drift and memory term near topography.
    1261                       term_1_2(n) = 0.0_wp
    1262                    ELSE IF ( num_gp > 1 )  THEN
    1263    
    1264                       d_sum = 0.0_wp
    1265 !
    1266 !--                   Evaluation of the distances between the gridpoints
    1267 !--                   contributing to the interpolated values, and the particle
    1268 !--                   location
    1269                       DO  agp = 1, num_gp
    1270                          d_gp_pl(agp) = ( particles(n)%x-location(agp,1) )**2  &
    1271                                       + ( particles(n)%y-location(agp,2) )**2  &
    1272                                       + ( zv(n)-location(agp,3) )**2
    1273                          d_sum        = d_sum + d_gp_pl(agp)
    1274                       ENDDO
    1275    
    1276 !
    1277 !--                  Finally the interpolation can be carried out
    1278                       e_int(n)     = 0.0_wp
    1279                       diss_int(n)  = 0.0_wp
    1280                       de_dx_int(n) = 0.0_wp 
    1281                       de_dy_int(n) = 0.0_wp 
    1282                       de_dz_int(n) = 0.0_wp 
    1283                       DO  agp = 1, num_gp
    1284                          e_int(n)     = e_int(n)     + ( d_sum - d_gp_pl(agp) ) * &
    1285                                                 ei(agp) / ( (num_gp-1) * d_sum )
    1286                          diss_int(n)  = diss_int(n)  + ( d_sum - d_gp_pl(agp) ) * &
    1287                                              dissi(agp) / ( (num_gp-1) * d_sum )
    1288                          de_dx_int(n) = de_dx_int(n) + ( d_sum - d_gp_pl(agp) ) * &
    1289                                             de_dxi(agp) / ( (num_gp-1) * d_sum )
    1290                          de_dy_int(n) = de_dy_int(n) + ( d_sum - d_gp_pl(agp) ) * &
    1291                                             de_dyi(agp) / ( (num_gp-1) * d_sum )
    1292                          de_dz_int(n) = de_dz_int(n) + ( d_sum - d_gp_pl(agp) ) * &
    1293                                             de_dzi(agp) / ( (num_gp-1) * d_sum )
    1294                       ENDDO
    1295    
    1296                    ENDIF
    1297                    e_int(n)     = MAX( 1E-20_wp, e_int(n)     )
    1298                    diss_int(n)  = MAX( 1E-20_wp, diss_int(n)  )
    1299                    de_dx_int(n) = MAX( 1E-20_wp, de_dx_int(n) )
    1300                    de_dy_int(n) = MAX( 1E-20_wp, de_dy_int(n) )
    1301                    de_dz_int(n) = MAX( 1E-20_wp, de_dz_int(n) )
    1302 !
    1303 !--                Set flag for stochastic equation which disables calculation
    1304 !--                of drift and memory term near topography.
    1305                    term_1_2(n) = 0.0_wp
    1306                 ENDIF
    1307              ENDDO
    1308           ENDDO
    1309        ENDIF
     685          ENDIF
     686       ENDDO
    1310687
    1311688       DO nb = 0,7
  • 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
  • palm/trunk/SOURCE/lpm_init.f90

    r2696 r2698  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Grid indices passed to lpm_boundary_conds. (responsible Philipp Thiele)
    2323!
    2424! Former revisions:
     
    236236
    237237    USE surface_mod,                                                           &
    238         ONLY:  get_topography_top_index, surf_def_h, surf_lsm_h, surf_usm_h
     238        ONLY:  get_topography_top_index_ji, surf_def_h, surf_lsm_h, surf_usm_h
    239239
    240240    IMPLICIT NONE
     
    757757!--                            Determine surface level. Therefore, check for
    758758!--                            upward-facing wall on w-grid.
    759                                k_surf = get_topography_top_index( jp, ip, 'w' )
     759                               k_surf = get_topography_top_index_ji( jp, ip, 'w' )
    760760
    761761                               IF ( seed_follows_topography )  THEN
     
    930930!--             Identify particles located outside the model domain and reflect
    931931!--             or absorb them if necessary.
    932                 CALL lpm_boundary_conds( 'bottom/top' )
     932                CALL lpm_boundary_conds( 'bottom/top', i, j, k )
    933933!
    934934!--             Furthermore, remove particles located in topography. Note, as
  • palm/trunk/SOURCE/microphysics_mod.f90

    r2696 r2698  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Bugfix in get_topography_top_index
    2323!
    2424! Former revisions:
     
    928928
    929929       USE surface_mod,                                                        &
    930            ONLY:  get_topography_top_index
     930           ONLY:  get_topography_top_index_ji
    931931
    932932
     
    945945!
    946946!--          Determine vertical index of topography top
    947              k_wall = get_topography_top_index( j, i, 's' )
     947             k_wall = get_topography_top_index_ji( j, i, 's' )
    948948             DO  k = nzb+1, nzt
    949949!
     
    23192319
    23202320       USE surface_mod,                                                        &
    2321            ONLY:  get_topography_top_index
     2321           ONLY:  get_topography_top_index_ji
    23222322
    23232323
     
    23342334!
    23352335!--    Determine vertical index of topography top
    2336        k_wall = get_topography_top_index( j, i, 's' )
     2336       k_wall = get_topography_top_index_ji( j, i, 's' )
    23372337       DO  k = nzb+1, nzt
    23382338!
  • palm/trunk/SOURCE/plant_canopy_model_mod.f90

    r2696 r2698  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Bugfix in get_topography_top_index
    2323!
    2424! Former revisions:
     
    151151
    152152    USE surface_mod,                                                           &
    153         ONLY:  get_topography_top_index
     153        ONLY:  get_topography_top_index_ji
    154154
    155155
     
    400400               DO  j = nys, nyn
    401401                  IF ( pch_index_ji(j,i) /= 0 )  THEN
    402                      k_topo = get_topography_top_index( j, i, 's' )
     402                     k_topo = get_topography_top_index_ji( j, i, 's' )
    403403                     DO  k = k_topo, k_topo + pch_index_ji(j,i)
    404404                        local_pf(i,j,k) = pc_heating_rate(k-k_topo,j,i)
     
    415415               DO  j = nys, nyn
    416416                  IF ( pch_index_ji(j,i) /= 0 )  THEN
    417                      k_topo = get_topography_top_index( j, i, 's' )
     417                     k_topo = get_topography_top_index_ji( j, i, 's' )
    418418                     DO  k = k_topo, k_topo + pch_index_ji(j,i)
    419419                        local_pf(i,j,k) = lad_s(k-k_topo,j,i)
     
    803803!--          Check whether topography and local vegetation on top exceed
    804804!--          height of the model domain.
    805              k = get_topography_top_index( j, i, 's' )
     805             k = get_topography_top_index_ji( j, i, 's' )
    806806             IF ( k + pch_index_ji(j,i) >= nzt + 1 )  THEN
    807807                message_string =  'Local vegetation height on top of ' //      &
     
    11101110!
    11111111!--                Determine topography-top index on u-grid
    1112                    k_wall = get_topography_top_index( j, i, 'u' )
     1112                   k_wall = get_topography_top_index_ji( j, i, 'u' )
    11131113                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
    11141114
     
    11741174!
    11751175!--                Determine topography-top index on v-grid
    1176                    k_wall = get_topography_top_index( j, i, 'v' )
     1176                   k_wall = get_topography_top_index_ji( j, i, 'v' )
    11771177
    11781178                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
     
    12391239!
    12401240!--                Determine topography-top index on w-grid
    1241                    k_wall = get_topography_top_index( j, i, 'w' )
     1241                   k_wall = get_topography_top_index_ji( j, i, 'w' )
    12421242
    12431243                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i) - 1
     
    12911291!
    12921292!--                Determine topography-top index on scalar-grid
    1293                    k_wall = get_topography_top_index( j, i, 's' )
     1293                   k_wall = get_topography_top_index_ji( j, i, 's' )
    12941294
    12951295                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
     
    13081308!
    13091309!--                Determine topography-top index on scalar-grid
    1310                    k_wall = get_topography_top_index( j, i, 's' )
     1310                   k_wall = get_topography_top_index_ji( j, i, 's' )
    13111311
    13121312                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
     
    13381338!
    13391339!--                Determine topography-top index on scalar-grid
    1340                    k_wall = get_topography_top_index( j, i, 's' )
     1340                   k_wall = get_topography_top_index_ji( j, i, 's' )
    13411341
    13421342                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
     
    13671367!
    13681368!--                Determine topography-top index on scalar-grid
    1369                    k_wall = get_topography_top_index( j, i, 's' )
     1369                   k_wall = get_topography_top_index_ji( j, i, 's' )
    13701370
    13711371                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
     
    14611461!
    14621462!--          Determine topography-top index on u-grid
    1463              k_wall = get_topography_top_index( j, i, 'u' )
     1463             k_wall = get_topography_top_index_ji( j, i, 'u' )
    14641464             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
    14651465
     
    15211521!
    15221522!--          Determine topography-top index on v-grid
    1523              k_wall = get_topography_top_index( j, i, 'v' )
     1523             k_wall = get_topography_top_index_ji( j, i, 'v' )
    15241524
    15251525             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
     
    15811581!
    15821582!--          Determine topography-top index on w-grid
    1583              k_wall = get_topography_top_index( j, i, 'w' )
     1583             k_wall = get_topography_top_index_ji( j, i, 'w' )
    15841584
    15851585             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i) - 1
     
    16281628!
    16291629!--          Determine topography-top index on scalar grid
    1630              k_wall = get_topography_top_index( j, i, 's' )
     1630             k_wall = get_topography_top_index_ji( j, i, 's' )
    16311631
    16321632             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
     
    16411641!
    16421642!--          Determine topography-top index on scalar grid
    1643              k_wall = get_topography_top_index( j, i, 's' )
     1643             k_wall = get_topography_top_index_ji( j, i, 's' )
    16441644
    16451645             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
     
    16671667!
    16681668!--          Determine topography-top index on scalar grid
    1669              k_wall = get_topography_top_index( j, i, 's' )
     1669             k_wall = get_topography_top_index_ji( j, i, 's' )
    16701670
    16711671             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
     
    16931693!
    16941694!--          Determine topography-top index on scalar grid
    1695              k_wall = get_topography_top_index( j, i, 's' )
     1695             k_wall = get_topography_top_index_ji( j, i, 's' )
    16961696
    16971697             DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
  • palm/trunk/SOURCE/pmc_interface_mod.f90

    r2696 r2698  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Bugfix in get_topography_top_index
    2323!
    2424! Former revisions:
     
    265265
    266266    USE surface_mod,                                                            &
    267         ONLY:  get_topography_top_index, surf_def_h, surf_lsm_h, surf_usm_h
     267        ONLY:  get_topography_top_index_ji, surf_def_h, surf_lsm_h, surf_usm_h
    268268
    269269    IMPLICIT NONE
     
    14101410!--             Determine largest topography index on scalar grid
    14111411                nzt_topo_nestbc_l = MAX( nzt_topo_nestbc_l,                    &
    1412                                          get_topography_top_index( j, i, 's' ) )
     1412                                         get_topography_top_index_ji( j, i, 's' ) )
    14131413!
    14141414!--             Determine largest topography index on u grid
    14151415                nzt_topo_nestbc_l = MAX( nzt_topo_nestbc_l,                    &
    1416                                          get_topography_top_index( j, i, 'u' ) )
     1416                                         get_topography_top_index_ji( j, i, 'u' ) )
    14171417!
    14181418!--             Determine largest topography index on v grid
    14191419                nzt_topo_nestbc_l = MAX( nzt_topo_nestbc_l,                    &
    1420                                          get_topography_top_index( j, i, 'v' ) )
     1420                                         get_topography_top_index_ji( j, i, 'v' ) )
    14211421!
    14221422!--             Determine largest topography index on w grid
    14231423                nzt_topo_nestbc_l = MAX( nzt_topo_nestbc_l,                    &
    1424                                          get_topography_top_index( j, i, 'w' ) )
     1424                                         get_topography_top_index_ji( j, i, 'w' ) )
    14251425             ENDDO
    14261426          ENDDO
     
    14361436!--             Determine largest topography index on scalar grid
    14371437                nzt_topo_nestbc_r = MAX( nzt_topo_nestbc_r,                    &
    1438                                          get_topography_top_index( j, i, 's' ) )
     1438                                         get_topography_top_index_ji( j, i, 's' ) )
    14391439!
    14401440!--             Determine largest topography index on u grid
    14411441                nzt_topo_nestbc_r = MAX( nzt_topo_nestbc_r,                    &
    1442                                          get_topography_top_index( j, i, 'u' ) )
     1442                                         get_topography_top_index_ji( j, i, 'u' ) )
    14431443!
    14441444!--             Determine largest topography index on v grid
    14451445                nzt_topo_nestbc_r = MAX( nzt_topo_nestbc_r,                    &
    1446                                          get_topography_top_index( j, i, 'v' ) )
     1446                                         get_topography_top_index_ji( j, i, 'v' ) )
    14471447!
    14481448!--             Determine largest topography index on w grid
    14491449                nzt_topo_nestbc_r = MAX( nzt_topo_nestbc_r,                    &
    1450                                          get_topography_top_index( j, i, 'w' ) )
     1450                                         get_topography_top_index_ji( j, i, 'w' ) )
    14511451          ENDDO
    14521452          nzt_topo_nestbc_r = nzt_topo_nestbc_r + 1
     
    14611461!--             Determine largest topography index on scalar grid
    14621462                nzt_topo_nestbc_s = MAX( nzt_topo_nestbc_s,                    &
    1463                                          get_topography_top_index( j, i, 's' ) )
     1463                                         get_topography_top_index_ji( j, i, 's' ) )
    14641464!
    14651465!--             Determine largest topography index on u grid
    14661466                nzt_topo_nestbc_s = MAX( nzt_topo_nestbc_s,                    &
    1467                                          get_topography_top_index( j, i, 'u' ) )
     1467                                         get_topography_top_index_ji( j, i, 'u' ) )
    14681468!
    14691469!--             Determine largest topography index on v grid
    14701470                nzt_topo_nestbc_s = MAX( nzt_topo_nestbc_s,                    &
    1471                                          get_topography_top_index( j, i, 'v' ) )
     1471                                         get_topography_top_index_ji( j, i, 'v' ) )
    14721472!
    14731473!--             Determine largest topography index on w grid
    14741474                nzt_topo_nestbc_s = MAX( nzt_topo_nestbc_s,                    &
    1475                                          get_topography_top_index( j, i, 'w' ) )
     1475                                         get_topography_top_index_ji( j, i, 'w' ) )
    14761476             ENDDO
    14771477          ENDDO
     
    14871487!--             Determine largest topography index on scalar grid
    14881488                nzt_topo_nestbc_n = MAX( nzt_topo_nestbc_n,                    &
    1489                                          get_topography_top_index( j, i, 's' ) )
     1489                                         get_topography_top_index_ji( j, i, 's' ) )
    14901490!
    14911491!--             Determine largest topography index on u grid
    14921492                nzt_topo_nestbc_n = MAX( nzt_topo_nestbc_n,                    &
    1493                                          get_topography_top_index( j, i, 'u' ) )
     1493                                         get_topography_top_index_ji( j, i, 'u' ) )
    14941494!
    14951495!--             Determine largest topography index on v grid
    14961496                nzt_topo_nestbc_n = MAX( nzt_topo_nestbc_n,                    &
    1497                                          get_topography_top_index( j, i, 'v' ) )
     1497                                         get_topography_top_index_ji( j, i, 'v' ) )
    14981498!
    14991499!--             Determine largest topography index on w grid
    15001500                nzt_topo_nestbc_n = MAX( nzt_topo_nestbc_n,                    &
    1501                                          get_topography_top_index( j, i, 'w' ) )
     1501                                         get_topography_top_index_ji( j, i, 'w' ) )
    15021502          ENDDO
    15031503          nzt_topo_nestbc_n = nzt_topo_nestbc_n + 1
     
    15521552!--          is part of the surfacetypes now. Set default roughness instead.
    15531553!--          Determine topography top index on u-grid
    1554              kb  = get_topography_top_index( j, i, 'u' )
     1554             kb  = get_topography_top_index_ji( j, i, 'u' )
    15551555             k   = kb + 1
    15561556             wall_index = kb
     
    15681568!
    15691569!--          Determine topography top index on v-grid
    1570              kb  = get_topography_top_index( j, i, 'v' )
     1570             kb  = get_topography_top_index_ji( j, i, 'v' )
    15711571             k   = kb + 1
    15721572             wall_index = kb
     
    16111611!--          to the present surface tpye.
    16121612!--          Determine topography top index on u-grid
    1613              kb  = get_topography_top_index( j, i, 'u' )
     1613             kb  = get_topography_top_index_ji( j, i, 'u' )
    16141614             k   = kb + 1
    16151615             wall_index = kb
     
    16271627!
    16281628!--          Determine topography top index on v-grid
    1629              kb  = get_topography_top_index( j, i, 'v' )
     1629             kb  = get_topography_top_index_ji( j, i, 'v' )
    16301630             k   = kb + 1
    16311631             wall_index = kb
     
    16671667!
    16681668!--          Determine topography top index on u-grid
    1669              kb  = get_topography_top_index( j, i, 'u' )
     1669             kb  = get_topography_top_index_ji( j, i, 'u' )
    16701670             k   = kb + 1
    16711671             wall_index = kb
     
    16831683!
    16841684!--          Determine topography top index on v-grid
    1685              kb  = get_topography_top_index( j, i, 'v' )
     1685             kb  = get_topography_top_index_ji( j, i, 'v' )
    16861686             k   = kb + 1
    16871687             wall_index = kb
     
    17231723!
    17241724!--          Determine topography top index on u-grid
    1725              kb  = get_topography_top_index( j, i, 'u' )
     1725             kb  = get_topography_top_index_ji( j, i, 'u' )
    17261726             k   = kb + 1
    17271727             wall_index = kb
     
    17391739!
    17401740!--          Determine topography top index on v-grid
    1741              kb  = get_topography_top_index( j, i, 'v' )
     1741             kb  = get_topography_top_index_ji( j, i, 'v' )
    17421742             k   = kb + 1
    17431743             wall_index = kb
     
    17791779!
    17801780!--             Determine lowest grid on outer grids for u and w.
    1781                 k_wall_u_ji   = get_topography_top_index( j,   0, 'u_out' )
    1782                 k_wall_u_ji_p = get_topography_top_index( j+1, 0, 'u_out' )
    1783                 k_wall_u_ji_m = get_topography_top_index( j-1, 0, 'u_out' )
    1784 
    1785                 k_wall_w_ji   = get_topography_top_index( j,   -1, 'w_out' )
    1786                 k_wall_w_ji_p = get_topography_top_index( j+1, -1, 'w_out' )
    1787                 k_wall_w_ji_m = get_topography_top_index( j-1, -1, 'w_out' )
     1781                k_wall_u_ji   = get_topography_top_index_ji( j,   0, 'u_out' )
     1782                k_wall_u_ji_p = get_topography_top_index_ji( j+1, 0, 'u_out' )
     1783                k_wall_u_ji_m = get_topography_top_index_ji( j-1, 0, 'u_out' )
     1784
     1785                k_wall_w_ji   = get_topography_top_index_ji( j,   -1, 'w_out' )
     1786                k_wall_w_ji_p = get_topography_top_index_ji( j+1, -1, 'w_out' )
     1787                k_wall_w_ji_m = get_topography_top_index_ji( j-1, -1, 'w_out' )
    17881788
    17891789                DO  k = nzb, nzt_topo_nestbc_l
     
    18691869!
    18701870!--             Determine lowest grid on outer grids for u and w.
    1871                 k_wall_u_ji   = get_topography_top_index( j,   i, 'u_out' )
    1872                 k_wall_u_ji_p = get_topography_top_index( j+1, i, 'u_out' )
    1873                 k_wall_u_ji_m = get_topography_top_index( j-1, i, 'u_out' )
    1874 
    1875                 k_wall_w_ji   = get_topography_top_index( j,   i, 'w_out' )
    1876                 k_wall_w_ji_p = get_topography_top_index( j+1, i, 'w_out' )
    1877                 k_wall_w_ji_m = get_topography_top_index( j-1, i, 'w_out' )
     1871                k_wall_u_ji   = get_topography_top_index_ji( j,   i, 'u_out' )
     1872                k_wall_u_ji_p = get_topography_top_index_ji( j+1, i, 'u_out' )
     1873                k_wall_u_ji_m = get_topography_top_index_ji( j-1, i, 'u_out' )
     1874
     1875                k_wall_w_ji   = get_topography_top_index_ji( j,   i, 'w_out' )
     1876                k_wall_w_ji_p = get_topography_top_index_ji( j+1, i, 'w_out' )
     1877                k_wall_w_ji_m = get_topography_top_index_ji( j-1, i, 'w_out' )
    18781878
    18791879                DO  k = nzb, nzt_topo_nestbc_r
     
    19521952!
    19531953!--             Determine lowest grid on outer grids for v and w.
    1954                 k_wall_v_ji   = get_topography_top_index( 0, i,   'v_out' )
    1955                 k_wall_v_ji_p = get_topography_top_index( 0, i+1, 'v_out' )
    1956                 k_wall_v_ji_m = get_topography_top_index( 0, i-1, 'v_out' )
    1957 
    1958                 k_wall_w_ji   = get_topography_top_index( -1, i,   'w_out' )
    1959                 k_wall_w_ji_p = get_topography_top_index( -1, i+1, 'w_out' )
    1960                 k_wall_w_ji_m = get_topography_top_index( -1, i-1, 'w_out' )
     1954                k_wall_v_ji   = get_topography_top_index_ji( 0, i,   'v_out' )
     1955                k_wall_v_ji_p = get_topography_top_index_ji( 0, i+1, 'v_out' )
     1956                k_wall_v_ji_m = get_topography_top_index_ji( 0, i-1, 'v_out' )
     1957
     1958                k_wall_w_ji   = get_topography_top_index_ji( -1, i,   'w_out' )
     1959                k_wall_w_ji_p = get_topography_top_index_ji( -1, i+1, 'w_out' )
     1960                k_wall_w_ji_m = get_topography_top_index_ji( -1, i-1, 'w_out' )
    19611961
    19621962                DO  k = nzb, nzt_topo_nestbc_s
     
    20412041!
    20422042!--             Determine lowest grid on outer grids for v and w.
    2043                 k_wall_v_ji   = get_topography_top_index( j, i,   'v_out' )
    2044                 k_wall_v_ji_p = get_topography_top_index( j, i+1, 'v_out' )
    2045                 k_wall_v_ji_m = get_topography_top_index( j, i-1, 'v_out' )
    2046 
    2047                 k_wall_w_ji   = get_topography_top_index( j, i,   'w_out' )
    2048                 k_wall_w_ji_p = get_topography_top_index( j, i+1, 'w_out' )
    2049                 k_wall_w_ji_m = get_topography_top_index( j, i-1, 'w_out' )
     2043                k_wall_v_ji   = get_topography_top_index_ji( j, i,   'v_out' )
     2044                k_wall_v_ji_p = get_topography_top_index_ji( j, i+1, 'v_out' )
     2045                k_wall_v_ji_m = get_topography_top_index_ji( j, i-1, 'v_out' )
     2046
     2047                k_wall_w_ji   = get_topography_top_index_ji( j, i,   'w_out' )
     2048                k_wall_w_ji_p = get_topography_top_index_ji( j, i+1, 'w_out' )
     2049                k_wall_w_ji_m = get_topography_top_index_ji( j, i-1, 'w_out' )
    20502050
    20512051                DO  k = nzb, nzt_topo_nestbc_n
     
    26432643          i = nxl - 1
    26442644          DO  j = nysg, nyng
    2645              k_wall = get_topography_top_index( j, i, 's' )
     2645             k_wall = get_topography_top_index_ji( j, i, 's' )
    26462646
    26472647             DO  k = k_wall + 1, nzt
     
    26642664          i = nxr + 1
    26652665          DO  j = nysg, nyng
    2666              k_wall = get_topography_top_index( j, i, 's' )
     2666             k_wall = get_topography_top_index_ji( j, i, 's' )
    26672667
    26682668             DO  k = k_wall + 1, nzt
     
    26852685          j = nys - 1
    26862686          DO  i = nxlg, nxrg
    2687              k_wall = get_topography_top_index( j, i, 's' )
     2687             k_wall = get_topography_top_index_ji( j, i, 's' )
    26882688             
    26892689             DO  k = k_wall + 1, nzt
     
    27062706          j = nyn + 1
    27072707          DO  i = nxlg, nxrg
    2708              k_wall = get_topography_top_index( j, i, 's' )
     2708             k_wall = get_topography_top_index_ji( j, i, 's' )
    27092709
    27102710             DO  k = k_wall + 1, nzt
     
    27292729!
    27302730!--          Determine vertical index for local topography top
    2731              k_wall = get_topography_top_index( j, i, 's' )
     2731             k_wall = get_topography_top_index_ji( j, i, 's' )
    27322732
    27332733             kc     = kco(k) + 1
     
    31963196!
    31973197!--             Determine vertical index of topography top at grid point (j,i)
    3198                 k_wall = get_topography_top_index( j, i, TRIM ( var ) )
     3198                k_wall = get_topography_top_index_ji( j, i, TRIM ( var ) )
    31993199!
    32003200!--             kbc is the first coarse-grid point above the surface
     
    32273227!
    32283228!--              Determine vertical index of topography top at grid point (j,i)
    3229                  k_wall = get_topography_top_index( j, i, 'w' )
     3229                 k_wall = get_topography_top_index_ji( j, i, 'w' )
    32303230
    32313231                 f(k_wall,j,i) = 0.0_wp
     
    43694369!
    43704370!--      Determine vertical index of topography top at grid point (j,i)
    4371          k_wall = get_topography_top_index( j, i, TRIM( var ) )
     4371         k_wall = get_topography_top_index_ji( j, i, TRIM( var ) )
    43724372
    43734373         DO  k = k_wall, nzt+1
     
    43954395!
    43964396!--         Determine vertical index of topography top at grid point (j,i)
    4397             k_wall = get_topography_top_index( j, i, TRIM ( var ) )
     4397            k_wall = get_topography_top_index_ji( j, i, TRIM ( var ) )
    43984398
    43994399            k = k_wall+1
     
    44204420!
    44214421!--            Determine vertical index of topography top at grid point (j,i)
    4422                k_wall = get_topography_top_index( j, i, TRIM ( var ) )
     4422               k_wall = get_topography_top_index_ji( j, i, TRIM ( var ) )
    44234423
    44244424               DO  k = k_wall+1, nzt_topo_nestbc
     
    44454445!
    44464446!--            Determine vertical index of topography top at grid point (j,i)
    4447                k_wall = get_topography_top_index( j, i, TRIM ( var ) )
     4447               k_wall = get_topography_top_index_ji( j, i, TRIM ( var ) )
    44484448
    44494449               k = k_wall + 1
     
    44764476!
    44774477!--            Determine vertical index of topography top at grid point (j,i)
    4478                k_wall = get_topography_top_index( j, i, 's' )
     4478               k_wall = get_topography_top_index_ji( j, i, 's' )
    44794479
    44804480               DO  k = k_wall, nzt + 1
     
    44864486!
    44874487!--            Determine vertical index of topography top at grid point (j,i)
    4488                k_wall = get_topography_top_index( j, i, 's' )
     4488               k_wall = get_topography_top_index_ji( j, i, 's' )
    44894489
    44904490               DO  k = k_wall, nzt+1
     
    45934593!
    45944594!--      Determine vertical index of topography top at grid point (j,i)
    4595          k_wall = get_topography_top_index( j, i, TRIM( var ) )
     4595         k_wall = get_topography_top_index_ji( j, i, TRIM( var ) )
    45964596
    45974597         DO  k = k_wall, nzt+1
     
    46194619!
    46204620!--         Determine vertical index of topography top at grid point (j,i)
    4621             k_wall = get_topography_top_index( j, i, TRIM( var ) )
     4621            k_wall = get_topography_top_index_ji( j, i, TRIM( var ) )
    46224622
    46234623            k = k_wall + 1
     
    46434643!
    46444644!--            Determine vertical index of topography top at grid point (j,i)
    4645                k_wall = get_topography_top_index( j, i, TRIM( var ) )
     4645               k_wall = get_topography_top_index_ji( j, i, TRIM( var ) )
    46464646
    46474647               DO  k = k_wall, nzt_topo_nestbc
     
    46704670!
    46714671!--            Determine vertical index of topography top at grid point (j,i)
    4672                k_wall = get_topography_top_index( j, i, TRIM( var ) )
     4672               k_wall = get_topography_top_index_ji( j, i, TRIM( var ) )
    46734673
    46744674               k = k_wall + 1
     
    47014701!
    47024702!--            Determine vertical index of topography top at grid point (j,i)
    4703                k_wall = get_topography_top_index( j, i, 's' )
     4703               k_wall = get_topography_top_index_ji( j, i, 's' )
    47044704               DO  k = k_wall, nzt+1
    47054705                  f(k,j,i) = tkefactor_s(k,i) * f(k,j,i)
     
    47104710!
    47114711!--            Determine vertical index of topography top at grid point (j,i)
    4712                k_wall = get_topography_top_index( j, i, 's' )
     4712               k_wall = get_topography_top_index_ji( j, i, 's' )
    47134713               DO  k = k_wall, nzt+1
    47144714                  f(k,j,i) = tkefactor_n(k,i) * f(k,j,i)
     
    48914891!
    48924892!--       Determine vertical index of topography top at grid point (j,i)
    4893           k_wall = get_topography_top_index( j, i, TRIM( var ) )
     4893          k_wall = get_topography_top_index_ji( j, i, TRIM( var ) )
    48944894
    48954895          DO  k = k_wall, nzt+1
     
    49704970!
    49714971!--       Determine vertical index of topography top at grid point (j,i)
    4972           k_wall = get_topography_top_index( j, i, TRIM( var ) )
     4972          k_wall = get_topography_top_index_ji( j, i, TRIM( var ) )
    49734973
    49744974          DO  k = k_wall, nzt+1
  • palm/trunk/SOURCE/radiation_model_mod.f90

    r2696 r2698  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Bugfix in get_topography_top_index
    2323!
    2424! Former revisions:
     
    272272
    273273    USE surface_mod,                                                           &
    274         ONLY:  get_topography_top_index, surf_def_h, surf_def_v, surf_lsm_h,   &
     274        ONLY:  get_topography_top_index, get_topography_top_index_ji,          &
     275               surf_def_h, surf_def_v, surf_lsm_h,                             &
    275276               surf_lsm_v, surf_type, surf_usm_h, surf_usm_v
    276277
     
    32443245!
    32453246!--             Obtain topography top index (lower bound of RRTMG)
    3246                 k_topo = get_topography_top_index( j, i, 's' )
     3247                k_topo = get_topography_top_index_ji( j, i, 's' )
    32473248
    32483249                IF ( lw_radiation )  THEN
     
    44114412                  pc_heating_rate, lad_s, prototype_lad, usm_lad_rma       
    44124413       
    4413        USE surface_mod,                                                        &
    4414            ONLY:  get_topography_top_index, surf_lsm_h, surf_lsm_v, surf_usm_h,&
    4415                   surf_usm_v
    4416 
    44174414       IMPLICIT NONE
    44184415
     
    44484445!
    44494446!--                Find topography top index
    4450                    k_topo = get_topography_top_index( j, i, 's' )
     4447                   k_topo = get_topography_top_index_ji( j, i, 's' )
    44514448
    44524449                   DO k = nzt+1, 0, -1
     
    45534550                DO  j = ijdb(3,ids), ijdb(4,ids)
    45544551
    4555                    k_topo  = get_topography_top_index( j, i, 's' )
    4556                    k_topo2 = get_topography_top_index( j-jdir(ids), i-idir(ids), 's' )
     4552                   k_topo  = get_topography_top_index_ji( j, i, 's' )
     4553                   k_topo2 = get_topography_top_index_ji( j-jdir(ids), i-idir(ids), 's' )
    45574554
    45584555
     
    45724569             DO j = nys, nyn
    45734570!--              Find topography top index
    4574                  k_topo = get_topography_top_index( j, i, 's' )
     4571                 k_topo = get_topography_top_index_ji( j, i, 's' )
    45754572                 k = nzut - k_topo
    45764573                 nsurfl = nsurfl + 6 * k
     
    45934590!
    45944591!--                Find topography top index
    4595                    k_topo = get_topography_top_index( j, i, 's' )
     4592                   k_topo = get_topography_top_index_ji( j, i, 's' )
    45964593
    45974594                   DO k = k_topo + 1, pct(j,i)
     
    47094706               DO i = ijdb(1,ids), ijdb(2,ids)
    47104707                   DO j = ijdb(3,ids), ijdb(4,ids)
    4711                        k_topo  = get_topography_top_index( j, i, 's' )
    4712                        k_topo2 = get_topography_top_index( j-jdir(ids), i-idir(ids), 's' )
     4708                       k_topo  = get_topography_top_index_ji( j, i, 's' )
     4709                       k_topo2 = get_topography_top_index_ji( j-jdir(ids), i-idir(ids), 's' )
    47134710
    47144711                       DO k = MAX(k_topo,k_topo2)+1, nzut
     
    47274724          DO i = nxl, nxr
    47284725             DO j = nys, nyn
    4729                 k_topo = get_topography_top_index( j, i, 's' )
     4726                k_topo = get_topography_top_index_ji( j, i, 's' )
    47304727
    47314728!--             add upward surface
     
    47464743          DO i = nxl, nxr
    47474744             DO j = nys, nyn
    4748                 k_topo = get_topography_top_index( j, i, 's' )
     4745                k_topo = get_topography_top_index_ji( j, i, 's' )
    47494746!--             north
    47504747                IF ( j /= ny ) THEN
     
    47524749                   jr = min(max(j-jdir(ids),0),ny)
    47534750                   ir = min(max(i-idir(ids),0),nx)
    4754                    k_topo2 = get_topography_top_index( jr, ir, 's' )
     4751                   k_topo2 = get_topography_top_index_ji( jr, ir, 's' )
    47554752                   DO k = MAX(k_topo,k_topo2)+1, nzut
    47564753                      isurf = isurf + 1
     
    47634760                   jr = min(max(j-jdir(ids),0),ny)
    47644761                   ir = min(max(i-idir(ids),0),nx)
    4765                    k_topo2 = get_topography_top_index( jr, ir, 's' )
     4762                   k_topo2 = get_topography_top_index_ji( jr, ir, 's' )
    47664763
    47674764                   DO k = MAX(k_topo,k_topo2)+1, nzut
     
    47754772                   jr = min(max(j-jdir(ids),0),ny)
    47764773                   ir = min(max(i-idir(ids),0),nx)
    4777                    k_topo2 = get_topography_top_index( jr, ir, 's' )
     4774                   k_topo2 = get_topography_top_index_ji( jr, ir, 's' )
    47784775
    47794776                   DO k = MAX(k_topo,k_topo2)+1, nzut
     
    47874784                   jr = min(max(j-jdir(ids),0),ny)
    47884785                   ir = min(max(i-idir(ids),0),nx)
    4789                    k_topo2 = get_topography_top_index( jr, ir, 's' )
     4786                   k_topo2 = get_topography_top_index_ji( jr, ir, 's' )
    47904787
    47914788                   DO k = MAX(k_topo,k_topo2)+1, nzut
     
    52555252!
    52565253!--             Following expression equals former kk = k - nzb_s_inner(j,i)
    5257                 kk = k - get_topography_top_index( j, i, 's' )  !- lad arrays are defined flat
     5254                kk = k - get_topography_top_index_ji( j, i, 's' )  !- lad arrays are defined flat
    52585255                pc_heating_rate(kk, j, i) = (pcbinsw(ipcgb) + pcbinlw(ipcgb)) &
    52595256                    * pchf_prep(k) * pt(k, j, i) !-- = dT/dt
     
    58565853            DO i = nxl, nxr
    58575854                DO j = nys, nyn
    5858                     k = get_topography_top_index( j, i, 's' )
     5855                    k = get_topography_top_index_ji( j, i, 's' )
    58595856
    58605857                    usm_lad(k:nzut, j, i) = lad_s(0:nzut-k, j, i)
  • palm/trunk/SOURCE/surface_mod.f90

    r2696 r2698  
    391391    PRIVATE
    392392
    393     INTERFACE get_topography_top_index
    394        MODULE PROCEDURE get_topography_top_index
    395        MODULE PROCEDURE get_topography_top_index_ji
    396     END  INTERFACE get_topography_top_index
    397 
    398393    INTERFACE init_bc
    399394       MODULE PROCEDURE init_bc
     
    431426!
    432427!-- Public subroutines and functions
    433     PUBLIC get_topography_top_index, init_bc, init_surfaces,                   &
     428    PUBLIC get_topography_top_index, get_topography_top_index_ji, init_bc,     &
     429           init_surfaces,                                                      &
    434430           init_surface_arrays, surface_read_restart_data,                     &
    435431           surface_restore_elements, surface_write_restart_data,               &
  • palm/trunk/SOURCE/turbulence_closure_mod.f90

    r2696 r2698  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Bugfix in get_topography_top_index
    2323!
    2424! Former revisions:
     
    870870
    871871    USE surface_mod,                                                           &
    872         ONLY:  get_topography_top_index
     872        ONLY:  get_topography_top_index_ji
    873873
    874874    IMPLICIT NONE
     
    972972          DO  i = nxlg, nxrg
    973973             DO  j = nysg, nyng
    974                 nz_s_shift = get_topography_top_index( j, i, 's' )
     974                nz_s_shift = get_topography_top_index_ji( j, i, 's' )
    975975
    976976                e(nz_s_shift:nzt+1,j,i)  =  e(0:nzt+1-nz_s_shift,j,i)
     
    991991          IF ( complex_terrain )  THEN
    992992             IF ( nxlg <= 0 .AND. nxrg >= 0 .AND. nysg <= 0 .AND. nyng >= 0 )  THEN
    993                 nz_s_shift_l = get_topography_top_index( 0, 0, 's' )
     993                nz_s_shift_l = get_topography_top_index_ji( 0, 0, 's' )
    994994             ELSE
    995995                nz_s_shift_l = 0
Note: See TracChangeset for help on using the changeset viewer.