Ignore:
Timestamp:
Mar 9, 2018 9:40:23 AM (3 years ago)
Author:
suehring
Message:

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

File:
1 edited

Legend:

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

    r2823 r2867  
    2525! -----------------
    2626! $Id$
     27! Revise mapping of 3D buildings onto onto orography.
     28!
     29! 2823 2018-02-20 15:31:45Z Giersch
    2730! Set boundary conditions for 3D topography in case of non-cyclic boundary
    2831! conditions
     
    835838    IMPLICIT NONE
    836839
    837     INTEGER(iwp) ::  i         !< running index along x-direction
    838     INTEGER(iwp) ::  j         !< running index along y-direction
    839     INTEGER(iwp) ::  k         !< running index along z-direction with respect to numeric grid
    840     INTEGER(iwp) ::  k2        !< running index along z-direction with respect to netcdf grid
    841     INTEGER(iwp) ::  k_surf    !< orography top index, used to map 3D buildings onto terrain
    842     INTEGER(iwp) ::  nr        !< index variable indication maximum terrain height for respective building ID
    843     INTEGER(iwp) ::  num_build !< counter for number of buildings
     840    INTEGER(iwp) ::  i                !< running index along x-direction
     841    INTEGER(iwp) ::  j                !< running index along y-direction
     842    INTEGER(iwp) ::  k                !< running index along z-direction with respect to numeric grid
     843    INTEGER(iwp) ::  k2               !< running index along z-direction with respect to netcdf grid
     844    INTEGER(iwp) ::  nr               !< index variable indication maximum terrain height for respective building ID
     845    INTEGER(iwp) ::  num_build        !< counter for number of buildings
     846    INTEGER(iwp) ::  topo_top_index   !< orography top index, used to map 3D buildings onto terrain
    844847
    845848    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  displace_dum        !< displacements of start addresses, used for MPI_ALLGATHERV
     
    992995!
    993996!--    Finally, determine maximumum terrain height occupied by the
    994 !--    respective building. First, on each PE locally.
     997!--    respective building.
    995998       ALLOCATE( oro_max_l(1:SIZE(build_ids_final)) )
    996999       ALLOCATE( oro_max(1:SIZE(build_ids_final))   )
     
    10131016#endif
    10141017!
    1015 !--    Map orography as well as buildings on grid.
     1018!--    Map orography as well as buildings onto grid.
    10161019!--    In case of ocean simulations, add an offset. 
    10171020       ocean_offset = MERGE( zw(0), 0.0_wp, ocean )
    10181021       DO  i = nxl, nxr
    10191022          DO  j = nys, nyn
     1023             topo_top_index = 0
    10201024             DO  k = nzb, nzt
    10211025!
     
    10301034                IF ( zu(k) - ocean_offset <= terrain_height_f%var(j,i) )  THEN
    10311035                    topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
    1032                     topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 1 )
     1036                    topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 1 )
     1037                    topo_top_index = topo_top_index + 1
    10331038                ENDIF
    10341039!
     
    10561061!
    10571062!--          Map 3D buildings onto terrain height. 
     1063!--          In case of any slopes, map building on top of maximum terrain
     1064!--          height covered by the building. In other words, extend
     1065!--          building down to the respective local terrain-surface height.
    10581066             IF ( buildings_f%lod == 2 )  THEN
    10591067                IF ( building_id_f%var(j,i) /= building_id_f%fill )  THEN
    10601068!
    1061 !--                Determine topo index occupied by terrain
    1062                    k_surf = MAXLOC( MERGE( 1, 0,                               &
    1063                                            .NOT. BTEST( topo_3d(:,j,i), 0 ) ), &
    1064                                     DIM = 1 ) - 1
    1065                    k2 = nzb+1
    1066                    DO k = k_surf + 1, nzt + 1
     1069!--                Determine index for maximum-terrain-height array.
     1070                   nr = MINLOC( ABS( build_ids_final -                         &
     1071                                     building_id_f%var(j,i) ), DIM = 1 )
     1072!
     1073!--                Extend building down to the terrain surface.
     1074                   k2 = topo_top_index
     1075                   DO k = topo_top_index + 1, nzt + 1     
     1076                      IF ( zu(k) - ocean_offset <= oro_max(nr) )  THEN
     1077                         topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
     1078                         topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
     1079                         k2             = k2 + 1
     1080                      ENDIF
     1081                   ENDDO   
     1082                   topo_top_index = k2       
     1083!
     1084!--                Now, map building on top.
     1085                   k2 = 0
     1086                   DO k = topo_top_index, nzt + 1
    10671087                      IF ( k2 <= buildings_f%nz-1 )  THEN
    10681088                         IF ( buildings_f%var_3d(k2,j,i) == 1 )  THEN
    10691089                            topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
     1090                            topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 1 )
    10701091                            topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
    10711092                         ENDIF
     
    20932114    ENDDO
    20942115
    2095     CALL exchange_horiz_int ( wall_flags_0, nys, nyn, nxl, nxr, nzt, nbgp )
     2116    CALL exchange_horiz_int( wall_flags_0, nys, nyn, nxl, nxr, nzt, nbgp )
    20962117!
    20972118!-- Set outer array for scalars to mask near-surface grid points in
     
    22972318       ENDDO
    22982319    ENDIF
    2299 
    23002320!
    23012321!-- Exchange ghost points for wall flags
Note: See TracChangeset for help on using the changeset viewer.