Changeset 2867 for palm

Timestamp:

Mar 9, 2018 9:40:23 AM (7 years ago)

Author:

suehring

Message:

Revise mapping of 3D-buildings onto orography; further bugfix concering call of user_init (see commit -r 2865)

Location:

palm/trunk/SOURCE

Files:

• init_3d_model.f90 (modified) (3 diffs)
• init_grid.f90 (modified) (8 diffs)

palm/trunk/SOURCE/init_3d_model.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Further bugfix concerning call of user_init. 
+! 
+! 2864 2018-03-08 11:57:45Z suehring
 ! Bugfix, move call of user_init in front of initialization of grid-point 
 ! arrays
@@ -2119 +2122 @@
     rmask(:,nxlg:nxl-1,:) = 0.0_wp;  rmask(:,nxr+1:nxrg,:) = 0.0_wp
     rmask(nysg:nys-1,:,:) = 0.0_wp;  rmask(nyn+1:nyng,:,:) = 0.0_wp
-
+!
+!-- User-defined initializing actions
+    CALL user_init
 !
 !-- To do: New concept for these non-topography grid points!
@@ -2177 +2182 @@
        ENDDO
     ENDDO
-
-!
-!-- User-defined initializing actions
-    CALL user_init
 !
 !-- Initialize arrays encompassing number of grid-points in inner and outer

palm/trunk/SOURCE/init_grid.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Revise mapping of 3D buildings onto onto orography.
+! 
+! 2823 2018-02-20 15:31:45Z Giersch
 ! Set boundary conditions for 3D topography in case of non-cyclic boundary 
 ! conditions
@@ -835 +838 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  i         !< running index along x-direction
-    INTEGER(iwp) ::  j         !< running index along y-direction
-    INTEGER(iwp) ::  k         !< running index along z-direction with respect to numeric grid
-    INTEGER(iwp) ::  k2        !< running index along z-direction with respect to netcdf grid
-    INTEGER(iwp) ::  k_surf    !< orography top index, used to map 3D buildings onto terrain
-    INTEGER(iwp) ::  nr        !< index variable indication maximum terrain height for respective building ID
-    INTEGER(iwp) ::  num_build !< counter for number of buildings
+    INTEGER(iwp) ::  i                !< running index along x-direction
+    INTEGER(iwp) ::  j                !< running index along y-direction
+    INTEGER(iwp) ::  k                !< running index along z-direction with respect to numeric grid
+    INTEGER(iwp) ::  k2               !< running index along z-direction with respect to netcdf grid
+    INTEGER(iwp) ::  nr               !< index variable indication maximum terrain height for respective building ID
+    INTEGER(iwp) ::  num_build        !< counter for number of buildings
+    INTEGER(iwp) ::  topo_top_index   !< orography top index, used to map 3D buildings onto terrain
 
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  displace_dum        !< displacements of start addresses, used for MPI_ALLGATHERV
@@ -992 +995 @@
 !
 !--    Finally, determine maximumum terrain height occupied by the 
-!--    respective building. First, on each PE locally. 
+!--    respective building. 
        ALLOCATE( oro_max_l(1:SIZE(build_ids_final)) )
        ALLOCATE( oro_max(1:SIZE(build_ids_final))   )
@@ -1013 +1016 @@
 #endif
 !
-!--    Map orography as well as buildings on grid. 
+!--    Map orography as well as buildings onto grid. 
 !--    In case of ocean simulations, add an offset.  
        ocean_offset = MERGE( zw(0), 0.0_wp, ocean )
        DO  i = nxl, nxr
           DO  j = nys, nyn
+             topo_top_index = 0
              DO  k = nzb, nzt
 !
@@ -1030 +1034 @@
                 IF ( zu(k) - ocean_offset <= terrain_height_f%var(j,i) )  THEN
                     topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
-                    topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 1 ) 
+                    topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 1 )
+                    topo_top_index = topo_top_index + 1
                 ENDIF
 !
@@ -1056 +1061 @@
 !
 !--          Map 3D buildings onto terrain height.  
+!--          In case of any slopes, map building on top of maximum terrain 
+!--          height covered by the building. In other words, extend 
+!--          building down to the respective local terrain-surface height. 
              IF ( buildings_f%lod == 2 )  THEN
                 IF ( building_id_f%var(j,i) /= building_id_f%fill )  THEN
 !
-!--                Determine topo index occupied by terrain
-                   k_surf = MAXLOC( MERGE( 1, 0,                               &
-                                           .NOT. BTEST( topo_3d(:,j,i), 0 ) ), &
-                                    DIM = 1 ) - 1
-                   k2 = nzb+1
-                   DO k = k_surf + 1, nzt + 1
+!--                Determine index for maximum-terrain-height array.
+                   nr = MINLOC( ABS( build_ids_final -                         &
+                                     building_id_f%var(j,i) ), DIM = 1 )
+!
+!--                Extend building down to the terrain surface.
+                   k2 = topo_top_index
+                   DO k = topo_top_index + 1, nzt + 1      
+                      IF ( zu(k) - ocean_offset <= oro_max(nr) )  THEN
+                         topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
+                         topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
+                         k2             = k2 + 1
+                      ENDIF
+                   ENDDO    
+                   topo_top_index = k2       
+!
+!--                Now, map building on top. 
+                   k2 = 0
+                   DO k = topo_top_index, nzt + 1
                       IF ( k2 <= buildings_f%nz-1 )  THEN
                          IF ( buildings_f%var_3d(k2,j,i) == 1 )  THEN
                             topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
+                            topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 1 )
                             topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
                          ENDIF
@@ -2093 +2114 @@
     ENDDO
 
-    CALL exchange_horiz_int ( wall_flags_0, nys, nyn, nxl, nxr, nzt, nbgp )
+    CALL exchange_horiz_int( wall_flags_0, nys, nyn, nxl, nxr, nzt, nbgp )
 !
 !-- Set outer array for scalars to mask near-surface grid points in 
@@ -2297 +2318 @@
        ENDDO
     ENDIF
-
 !
 !-- Exchange ghost points for wall flags