Changeset 3051 for palm/trunk/SOURCE/init_grid.f90

Timestamp:

May 30, 2018 5:43:55 PM (6 years ago)

Author:

suehring

Message:

Speed-up NetCDF input; Revise NetCDF-input routines and remove input via io-blocks; Temporarily revoke renaming of input variables in dynamic driver; More detailed error messages created; Bugfix in mapping 3D buildings; Bugfix in land-surface model at pavement surfaces; Bugfix in initialization with inifor

File:

• palm/trunk/SOURCE/init_grid.f90 (modified) (6 diffs)

palm/trunk/SOURCE/init_grid.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Minor bugfix concerning mapping 3D buildings on top of terrain
+! 
+! 3045 2018-05-28 07:55:41Z Giersch
 ! Error messages revised
 ! 
@@ -926 +929 @@
     topo_3d(nzb,:,:) = IBCLR( topo_3d(nzb,:,:), 0 )
 !
+!-- In order to map topography on PALM grid also in case of ocean simulations,
+!-- pre-calculate an offset value.
+    ocean_offset = MERGE( zw(0), 0.0_wp, ocean )
+!
 !-- Reference buildings on top of orography. This is not necessary
 !-- if topography is read from ASCII file as no distinction between buildings
@@ -1038 +1045 @@
 
 !
-!--       Finally, determine maximumum terrain height occupied by the 
-!--       respective building. 
+!--       Determine maximumum terrain height occupied by the respective 
+!--       building and temporalily store on oro_max 
           ALLOCATE( oro_max_l(1:SIZE(build_ids_final)) )
           ALLOCATE( oro_max(1:SIZE(build_ids_final))   )
@@ -1059 +1066 @@
           oro_max = oro_max_l
 #endif
+!
+!--       Finally, determine discrete grid height of maximum orography occupied
+!--       by a building. Use all-or-nothing approach, i.e. a grid box is either
+          oro_max_l = 0.0
+          DO  nr = 1, SIZE(build_ids_final)
+             DO  k = nzb, nzt
+                IF ( zu(k) - ocean_offset <= oro_max(nr) )                     &
+                   oro_max_l = zw(k) - ocean_offset
+             ENDDO
+             oro_max = oro_max_l
+          ENDDO
        ENDIF
 !
 !--    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
@@ -1080 +1096 @@
                     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_top_index = topo_top_index + 1
+                    topo_top_index = k ! topo_top_index + 1
                 ENDIF
 !
@@ -1116 +1132 @@
                                      building_id_f%var(j,i) ), DIM = 1 )
 !
-!--                Extend building down to the terrain surface.
-                   k2 = topo_top_index
+!--                Extend building down to the terrain surface, i.e. fill-up
+!--                surface irregularities below a building. Note, oro_max
+!--                is already a discrete height according to the all-or-nothing
+!--                approach, i.e. grid box is either topography or atmosphere, 
+!--                terrain top is defined at upper bound of the grid box.
+!--                Hence, check for zw in this case. 
                    DO k = topo_top_index + 1, nzt + 1      
-                      IF ( zu(k) - ocean_offset <= oro_max(nr) )  THEN
+                      IF ( zw(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. 
+                   ENDDO        
+!
+!--                After surface irregularities are smoothen, determine lower
+!--                start index where building starts. 
+                   DO  k = nzb, nzt
+                      IF ( zw(k) - ocean_offset <= oro_max(nr) )               &
+                         topo_top_index = k
+                   ENDDO
+!
+!--                Finally, map building on top.
                    k2 = 0
                    DO k = topo_top_index, nzt + 1