Ignore:
Timestamp:
Feb 23, 2009 1:03:18 PM (16 years ago)
Author:
raasch
Message:

further additions for clipping - still incomplete

File:
1 edited

Legend:

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

    r237 r242  
    44! Actual revisions:
    55! -----------------
     6! Clipping implemented.
    67! Polygon reduction for building and ground plate isosurface. Reduction level
    78! for buildings can be chosen with parameter cluster_size.
     
    6061    CHARACTER (LEN=80) ::  dvrp_file_local
    6162    INTEGER ::  cluster_mode, cluster_size_x, cluster_size_y, cluster_size_z, &
    62                 gradient_normals, i, j, k, l, m, pn, tv, vn
     63                gradient_normals, i, j, k, l, m, nx_dvrp_l, nx_dvrp_r,        &
     64                ny_dvrp_n, ny_dvrp_s, pn, tv, vn
    6365    LOGICAL ::  allocated
    6466    REAL(4) ::  center(3), cluster_alpha, distance, tmp_b, tmp_g, tmp_r, &
     
    7173                            dvrp_file_local_c,dvrp_host_c, &
    7274                            dvrp_password_c, dvrp_username_c, name_c
     75
     76!
     77!-- Set clipping to default (total domain), if not set by user
     78    IF ( clip_dvrp_l == 9999999.9 )  clip_dvrp_l = 0.0
     79    IF ( clip_dvrp_r == 9999999.9 )  clip_dvrp_r = ( nx + 1 ) * dx
     80    IF ( clip_dvrp_s == 9999999.9 )  clip_dvrp_s = 0.0
     81    IF ( clip_dvrp_n == 9999999.9 )  clip_dvrp_n = ( ny + 1 ) * dy
     82
     83!
     84!-- Calculate the clipping index limits
     85    nx_dvrp_l = clip_dvrp_l / dx
     86    nx_dvrp_r = clip_dvrp_r / dx
     87    ny_dvrp_s = clip_dvrp_s / dy
     88    ny_dvrp_n = clip_dvrp_n / dy
     89
     90    IF ( nx_dvrp_l < nxr  .AND.  nx_dvrp_r > nxl  .AND. &
     91         ny_dvrp_s < nyn  .AND.  ny_dvrp_n > nys )  THEN
     92
     93       dvrp_overlap = .TRUE.
     94       nxl_dvrp = MAX( nxl, nx_dvrp_l )
     95       nxr_dvrp = MIN( nxr, nx_dvrp_r )
     96       nys_dvrp = MAX( nys, ny_dvrp_s )
     97       nyn_dvrp = MIN( nyn, ny_dvrp_n )
     98
     99       IF ( nxl_dvrp == nxl  .AND.  nxr_dvrp == nxr  .AND.  &
     100            nys_dvrp == nys  .AND.  nyn_dvrp == nyn )  THEN
     101          dvrp_total_overlap = .TRUE.
     102       ELSE
     103          dvrp_total_overlap = .FALSE.
     104       ENDIF
     105
     106    ELSE
     107!
     108!--    This subdomain does not overlap with the clipping area. Define an
     109!--    arbitrary (small) domain within in the clipping area.
     110       dvrp_overlap       = .FALSE.
     111       dvrp_total_overlap = .FALSE.
     112       nxl_dvrp = nx_dvrp_l
     113       nxr_dvrp = nxl_dvrp + 4
     114       nys_dvrp = ny_dvrp_s
     115       nyn_dvrp = nys_dvrp + 4
     116
     117    ENDIF
    73118
    74119!
     
    151196!
    152197!--       Compute center of domain and distance of camera from center
    153           center(1) = ( nx + 1.0 ) * dx * 0.5 * superelevation_x
    154           center(2) = ( ny + 1.0 ) * dy * 0.5 * superelevation_y
     198          center(1) = ( clip_dvrp_l + clip_dvrp_r ) * 0.5 * superelevation_x
     199          center(2) = ( clip_dvrp_s + clip_dvrp_n ) * 0.5 * superelevation_y
    155200          center(3) = ( zu(nz_do3d) - zu(nzb) ) * 0.5 * superelevation
    156           distance  = 1.5 * MAX( ( nx + 1.0 ) * dx * superelevation_x, &
    157                                  ( ny + 1.0 ) * dy * superelevation_y, &
     201          distance  = 1.5 * MAX( (clip_dvrp_r-clip_dvrp_l) * superelevation_x, &
     202                                 (clip_dvrp_n-clip_dvrp_s) * superelevation_y, &
    158203                                 ( zu(nz_do3d) - zu(nzb) ) * superelevation )
    159204
     
    200245          CALL DVRP_MATERIAL_RGB( m-1, 1, tmp_r, tmp_g, tmp_b, tmp_t )
    201246
    202           tmp_1 = 0.01;  tmp_2 = 0.0;  tmp_3 = 0.0;  tmp_4 = 0.0
    203           tmp_5 = (nx+1) * dx * superelevation_x
    204           tmp_6 = (ny+1) * dy * superelevation_y
     247          tmp_1 = 0.01;
     248          tmp_2 = clip_dvrp_l * superelevation_x
     249          tmp_3 = clip_dvrp_s * superelevation_y
     250          tmp_4 = 0.0
     251          tmp_5 = clip_dvrp_r * superelevation_x
     252          tmp_6 = clip_dvrp_n * superelevation_y
    205253          tmp_7 = zu(nz_do3d) * superelevation
    206254          CALL DVRP_BOUNDINGBOX( m-1, 1, tmp_1, tmp_2, tmp_3, tmp_4, tmp_5, &
     
    252300!--          Determine local gridpoint coordinates
    253301             IF ( .NOT. allocated )  THEN
    254                 ALLOCATE( xcoor_dvrp(nxl:nxr+1), ycoor_dvrp(nys:nyn+1), &
     302                ALLOCATE( xcoor_dvrp(nxl_dvrp:nxr_dvrp+1), &
     303                          ycoor_dvrp(nys_dvrp:nyn_dvrp+1), &
    255304                          zcoor_dvrp(nzb:nz_do3d) )
    256305                allocated = .TRUE.
    257306
    258                 DO  i = nxl, nxr+1
     307                DO  i = nxl_dvrp, nxr_dvrp+1
    259308                   xcoor_dvrp(i) = i * dx * superelevation_x
    260309                ENDDO
    261                 DO  j = nys, nyn+1
     310                DO  j = nys_dvrp, nyn_dvrp+1
    262311                   ycoor_dvrp(j) = j * dy * superelevation_y
    263312                ENDDO
    264313                zcoor_dvrp = zu(nzb:nz_do3d) * superelevation
    265                 nx_dvrp    = nxr+1 - nxl + 1
    266                 ny_dvrp    = nyn+1 - nys + 1
     314                nx_dvrp    = nxr_dvrp+1 - nxl_dvrp + 1
     315                ny_dvrp    = nyn_dvrp+1 - nys_dvrp + 1
    267316                nz_dvrp    = nz_do3d - nzb + 1
    268317             ENDIF
     
    279328!
    280329!--          Compute and plot isosurface in dvr-format
    281              ALLOCATE( local_pf(nxl:nxr+1,nys:nyn+1,nzb:nz_do3d) )
     330             ALLOCATE( local_pf(nxl_dvrp:nxr_dvrp+1,nys_dvrp:nyn_dvrp+1, &
     331                                nzb:nz_do3d) )
    282332             local_pf = 0.0
    283              DO  i = nxl, nxr+1
    284                 DO  j = nys, nyn+1
    285                    IF ( nzb_s_inner(j,i) > 0 )  THEN
     333             IF ( dvrp_overlap )  THEN
     334                DO  i = nxl_dvrp, nxr_dvrp+1
     335                   DO  j = nys_dvrp, nyn_dvrp+1
     336                      IF ( nzb_s_inner(j,i) > 0 )  THEN
    286337                         local_pf(i,j,nzb:nzb_s_inner(j,i)) = 1.0
    287338                      ENDIF
     339                   ENDDO
    288340                ENDDO
    289              ENDDO
     341             ENDIF
    290342
    291343             CALL DVRP_DATA( m-1, local_pf, 1, nx_dvrp, ny_dvrp, nz_dvrp, &
     
    369421!--       Determine local gridpoint coordinates
    370422          IF ( .NOT. allocated )  THEN
    371              ALLOCATE( xcoor_dvrp(nxl:nxr+1), ycoor_dvrp(nys:nyn+1), &
     423             ALLOCATE( xcoor_dvrp(nxl_dvrp:nxr_dvrp+1), &
     424                       ycoor_dvrp(nys_dvrp:nyn_dvrp+1), &
    372425                       zcoor_dvrp(nzb:nz_do3d) )
    373426             allocated = .TRUE.
    374427
    375              DO  i = nxl, nxr+1
     428             DO  i = nxl_dvrp, nxr_dvrp+1
    376429                xcoor_dvrp(i) = i * dx * superelevation_x
    377430             ENDDO
    378              DO  j = nys, nyn+1
     431             DO  j = nys_dvrp, nyn_dvrp+1
    379432                ycoor_dvrp(j) = j * dy * superelevation_y
    380433             ENDDO
    381434             zcoor_dvrp = zu(nzb:nz_do3d) * superelevation
    382              nx_dvrp    = nxr+1 - nxl + 1
    383              ny_dvrp    = nyn+1 - nys + 1
     435             nx_dvrp    = nxr_dvrp+1 - nxl_dvrp + 1
     436             ny_dvrp    = nyn_dvrp+1 - nys_dvrp + 1
    384437             nz_dvrp    = nz_do3d - nzb + 1
    385438          ENDIF
     
    396449!
    397450!--       Compute and plot isosurface in dvr-format
    398           ALLOCATE( local_pf(nxl:nxr+1,nys:nyn+1,nzb:nz_do3d) )
     451          ALLOCATE( local_pf(nxl_dvrp:nxr_dvrp+1,nys_dvrp:nyn_dvrp+1, &
     452                             nzb:nz_do3d) )
    399453          local_pf = 0.0
    400           local_pf(:,:,0) = 1.0
     454          IF (dvrp_overlap )  local_pf(:,:,0) = 1.0
    401455
    402456          CALL DVRP_DATA( m-1, local_pf, 1, nx_dvrp, ny_dvrp, nz_dvrp, &
     
    406460
    407461!
    408 !--       Always reduce the number of polygones
     462!--       Always reduce the number of polygones as much as possible
    409463          cluster_size_x = 5
    410464          cluster_size_y = 5
     
    541595!--    Determine local gridpoint coordinates
    542596       IF ( .NOT. allocated )  THEN
    543           ALLOCATE( xcoor_dvrp(nxl:nxr+1), ycoor_dvrp(nys:nyn+1), &
     597          ALLOCATE( xcoor_dvrp(nxl_dvrp:nxr_dvrp+1), &
     598                    ycoor_dvrp(nys_dvrp:nyn_dvrp+1), &
    544599                    zcoor_dvrp(nzb:nz_do3d) )
    545600          allocated = .TRUE.
    546601
    547           DO  i = nxl, nxr+1
     602          DO  i = nxl_dvrp, nxr_dvrp+1
    548603             xcoor_dvrp(i) = i * dx * superelevation_x
    549604          ENDDO
    550           DO  j = nys, nyn+1
     605          DO  j = nys_dvrp, nyn_dvrp+1
    551606             ycoor_dvrp(j) = j * dy * superelevation_y
    552607          ENDDO
    553608          zcoor_dvrp = zu(nzb:nz_do3d) * superelevation
    554           nx_dvrp    = nxr+1 - nxl + 1
    555           ny_dvrp    = nyn+1 - nys + 1
     609          nx_dvrp    = nxr_dvrp+1 - nxl_dvrp + 1
     610          ny_dvrp    = nyn_dvrp+1 - nys_dvrp + 1
    556611          nz_dvrp    = nz_do3d - nzb + 1
    557612       ENDIF
Note: See TracChangeset for help on using the changeset viewer.