Changeset 242 for palm/trunk/SOURCE/init_dvrp.f90

Timestamp:

Feb 23, 2009 1:03:18 PM (15 years ago)

Author:

raasch

Message:

further additions for clipping - still incomplete

File:

• palm/trunk/SOURCE/init_dvrp.f90 (modified) (11 diffs)

palm/trunk/SOURCE/init_dvrp.f90

@@ -4 +4 @@
 ! Actual revisions:
 ! -----------------
+! Clipping implemented.
 ! Polygon reduction for building and ground plate isosurface. Reduction level
 ! for buildings can be chosen with parameter cluster_size.
@@ -60 +61 @@
     CHARACTER (LEN=80) ::  dvrp_file_local
     INTEGER ::  cluster_mode, cluster_size_x, cluster_size_y, cluster_size_z, &
-                gradient_normals, i, j, k, l, m, pn, tv, vn
+                gradient_normals, i, j, k, l, m, nx_dvrp_l, nx_dvrp_r,        &
+                ny_dvrp_n, ny_dvrp_s, pn, tv, vn
     LOGICAL ::  allocated
     REAL(4) ::  center(3), cluster_alpha, distance, tmp_b, tmp_g, tmp_r, &
@@ -71 +73 @@
                             dvrp_file_local_c,dvrp_host_c, &
                             dvrp_password_c, dvrp_username_c, name_c
+
+!
+!-- Set clipping to default (total domain), if not set by user
+    IF ( clip_dvrp_l == 9999999.9 )  clip_dvrp_l = 0.0
+    IF ( clip_dvrp_r == 9999999.9 )  clip_dvrp_r = ( nx + 1 ) * dx
+    IF ( clip_dvrp_s == 9999999.9 )  clip_dvrp_s = 0.0
+    IF ( clip_dvrp_n == 9999999.9 )  clip_dvrp_n = ( ny + 1 ) * dy
+
+!
+!-- Calculate the clipping index limits
+    nx_dvrp_l = clip_dvrp_l / dx
+    nx_dvrp_r = clip_dvrp_r / dx
+    ny_dvrp_s = clip_dvrp_s / dy
+    ny_dvrp_n = clip_dvrp_n / dy
+
+    IF ( nx_dvrp_l < nxr  .AND.  nx_dvrp_r > nxl  .AND. &
+         ny_dvrp_s < nyn  .AND.  ny_dvrp_n > nys )  THEN
+
+       dvrp_overlap = .TRUE.
+       nxl_dvrp = MAX( nxl, nx_dvrp_l )
+       nxr_dvrp = MIN( nxr, nx_dvrp_r )
+       nys_dvrp = MAX( nys, ny_dvrp_s )
+       nyn_dvrp = MIN( nyn, ny_dvrp_n )
+
+       IF ( nxl_dvrp == nxl  .AND.  nxr_dvrp == nxr  .AND.  &
+            nys_dvrp == nys  .AND.  nyn_dvrp == nyn )  THEN
+          dvrp_total_overlap = .TRUE.
+       ELSE
+          dvrp_total_overlap = .FALSE.
+       ENDIF
+
+    ELSE
+!
+!--    This subdomain does not overlap with the clipping area. Define an
+!--    arbitrary (small) domain within in the clipping area.
+       dvrp_overlap       = .FALSE.
+       dvrp_total_overlap = .FALSE.
+       nxl_dvrp = nx_dvrp_l
+       nxr_dvrp = nxl_dvrp + 4
+       nys_dvrp = ny_dvrp_s
+       nyn_dvrp = nys_dvrp + 4
+
+    ENDIF
 
 !
@@ -151 +196 @@
 !
 !--       Compute center of domain and distance of camera from center
-          center(1) = ( nx + 1.0 ) * dx * 0.5 * superelevation_x
-          center(2) = ( ny + 1.0 ) * dy * 0.5 * superelevation_y
+          center(1) = ( clip_dvrp_l + clip_dvrp_r ) * 0.5 * superelevation_x
+          center(2) = ( clip_dvrp_s + clip_dvrp_n ) * 0.5 * superelevation_y
           center(3) = ( zu(nz_do3d) - zu(nzb) ) * 0.5 * superelevation
-          distance  = 1.5 * MAX( ( nx + 1.0 ) * dx * superelevation_x, &
-                                 ( ny + 1.0 ) * dy * superelevation_y, &
+          distance  = 1.5 * MAX( (clip_dvrp_r-clip_dvrp_l) * superelevation_x, &
+                                 (clip_dvrp_n-clip_dvrp_s) * superelevation_y, &
                                  ( zu(nz_do3d) - zu(nzb) ) * superelevation )
 
@@ -200 +245 @@
           CALL DVRP_MATERIAL_RGB( m-1, 1, tmp_r, tmp_g, tmp_b, tmp_t )
 
-          tmp_1 = 0.01;  tmp_2 = 0.0;  tmp_3 = 0.0;  tmp_4 = 0.0
-          tmp_5 = (nx+1) * dx * superelevation_x
-          tmp_6 = (ny+1) * dy * superelevation_y
+          tmp_1 = 0.01;
+          tmp_2 = clip_dvrp_l * superelevation_x
+          tmp_3 = clip_dvrp_s * superelevation_y
+          tmp_4 = 0.0
+          tmp_5 = clip_dvrp_r * superelevation_x
+          tmp_6 = clip_dvrp_n * superelevation_y
           tmp_7 = zu(nz_do3d) * superelevation
           CALL DVRP_BOUNDINGBOX( m-1, 1, tmp_1, tmp_2, tmp_3, tmp_4, tmp_5, &
@@ -252 +300 @@
 !--          Determine local gridpoint coordinates
              IF ( .NOT. allocated )  THEN
-                ALLOCATE( xcoor_dvrp(nxl:nxr+1), ycoor_dvrp(nys:nyn+1), &
+                ALLOCATE( xcoor_dvrp(nxl_dvrp:nxr_dvrp+1), &
+                          ycoor_dvrp(nys_dvrp:nyn_dvrp+1), &
                           zcoor_dvrp(nzb:nz_do3d) )
                 allocated = .TRUE.
 
-                DO  i = nxl, nxr+1
+                DO  i = nxl_dvrp, nxr_dvrp+1
                    xcoor_dvrp(i) = i * dx * superelevation_x
                 ENDDO
-                DO  j = nys, nyn+1
+                DO  j = nys_dvrp, nyn_dvrp+1
                    ycoor_dvrp(j) = j * dy * superelevation_y
                 ENDDO
                 zcoor_dvrp = zu(nzb:nz_do3d) * superelevation
-                nx_dvrp    = nxr+1 - nxl + 1
-                ny_dvrp    = nyn+1 - nys + 1
+                nx_dvrp    = nxr_dvrp+1 - nxl_dvrp + 1
+                ny_dvrp    = nyn_dvrp+1 - nys_dvrp + 1
                 nz_dvrp    = nz_do3d - nzb + 1
              ENDIF
@@ -279 +328 @@
 !
 !--          Compute and plot isosurface in dvr-format
-             ALLOCATE( local_pf(nxl:nxr+1,nys:nyn+1,nzb:nz_do3d) )
+             ALLOCATE( local_pf(nxl_dvrp:nxr_dvrp+1,nys_dvrp:nyn_dvrp+1, &
+                                nzb:nz_do3d) )
              local_pf = 0.0
-             DO  i = nxl, nxr+1
-                DO  j = nys, nyn+1
-                   IF ( nzb_s_inner(j,i) > 0 )  THEN
+             IF ( dvrp_overlap )  THEN
+                DO  i = nxl_dvrp, nxr_dvrp+1
+                   DO  j = nys_dvrp, nyn_dvrp+1
+                      IF ( nzb_s_inner(j,i) > 0 )  THEN
                          local_pf(i,j,nzb:nzb_s_inner(j,i)) = 1.0
                       ENDIF
+                   ENDDO
                 ENDDO
-             ENDDO
+             ENDIF
 
              CALL DVRP_DATA( m-1, local_pf, 1, nx_dvrp, ny_dvrp, nz_dvrp, &
@@ -369 +421 @@
 !--       Determine local gridpoint coordinates
           IF ( .NOT. allocated )  THEN
-             ALLOCATE( xcoor_dvrp(nxl:nxr+1), ycoor_dvrp(nys:nyn+1), &
+             ALLOCATE( xcoor_dvrp(nxl_dvrp:nxr_dvrp+1), &
+                       ycoor_dvrp(nys_dvrp:nyn_dvrp+1), &
                        zcoor_dvrp(nzb:nz_do3d) )
              allocated = .TRUE.
 
-             DO  i = nxl, nxr+1
+             DO  i = nxl_dvrp, nxr_dvrp+1
                 xcoor_dvrp(i) = i * dx * superelevation_x
              ENDDO
-             DO  j = nys, nyn+1
+             DO  j = nys_dvrp, nyn_dvrp+1
                 ycoor_dvrp(j) = j * dy * superelevation_y
              ENDDO
              zcoor_dvrp = zu(nzb:nz_do3d) * superelevation
-             nx_dvrp    = nxr+1 - nxl + 1
-             ny_dvrp    = nyn+1 - nys + 1
+             nx_dvrp    = nxr_dvrp+1 - nxl_dvrp + 1
+             ny_dvrp    = nyn_dvrp+1 - nys_dvrp + 1
              nz_dvrp    = nz_do3d - nzb + 1
           ENDIF
@@ -396 +449 @@
 !
 !--       Compute and plot isosurface in dvr-format
-          ALLOCATE( local_pf(nxl:nxr+1,nys:nyn+1,nzb:nz_do3d) )
+          ALLOCATE( local_pf(nxl_dvrp:nxr_dvrp+1,nys_dvrp:nyn_dvrp+1, &
+                             nzb:nz_do3d) )
           local_pf = 0.0
-          local_pf(:,:,0) = 1.0
+          IF (dvrp_overlap )  local_pf(:,:,0) = 1.0
 
           CALL DVRP_DATA( m-1, local_pf, 1, nx_dvrp, ny_dvrp, nz_dvrp, &
@@ -406 +460 @@
 
 !
-!--       Always reduce the number of polygones
+!--       Always reduce the number of polygones as much as possible
           cluster_size_x = 5
           cluster_size_y = 5
@@ -541 +595 @@
 !--    Determine local gridpoint coordinates
        IF ( .NOT. allocated )  THEN
-          ALLOCATE( xcoor_dvrp(nxl:nxr+1), ycoor_dvrp(nys:nyn+1), &
+          ALLOCATE( xcoor_dvrp(nxl_dvrp:nxr_dvrp+1), &
+                    ycoor_dvrp(nys_dvrp:nyn_dvrp+1), &
                     zcoor_dvrp(nzb:nz_do3d) )
           allocated = .TRUE.
 
-          DO  i = nxl, nxr+1
+          DO  i = nxl_dvrp, nxr_dvrp+1
              xcoor_dvrp(i) = i * dx * superelevation_x
           ENDDO
-          DO  j = nys, nyn+1
+          DO  j = nys_dvrp, nyn_dvrp+1
              ycoor_dvrp(j) = j * dy * superelevation_y
           ENDDO
           zcoor_dvrp = zu(nzb:nz_do3d) * superelevation
-          nx_dvrp    = nxr+1 - nxl + 1
-          ny_dvrp    = nyn+1 - nys + 1
+          nx_dvrp    = nxr_dvrp+1 - nxl_dvrp + 1
+          ny_dvrp    = nyn_dvrp+1 - nys_dvrp + 1
           nz_dvrp    = nz_do3d - nzb + 1
        ENDIF