Changeset 4189 for palm/trunk

Timestamp:

Aug 26, 2019 4:19:38 PM (5 years ago)

Author:

suehring

Message:

init_grid: Add check for proper setting of namelist parameter topography; Set bit flag to indicate land surfaces in case no topography is provided

File:

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

palm/trunk/SOURCE/init_grid.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! - Add check for proper setting of namelist parameter topography
+! - Set flag to indicate land surfaces in case no topography is provided
+! 
+! 4182 2019-08-22 15:20:23Z scharf
 ! Corrected "Former revisions" section
 ! 
@@ -1333 +1337 @@
 !--             Hence, in this case, also a building flag.
                 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_top_index = k ! topo_top_index + 1
+                   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 = k ! topo_top_index + 1
                 ENDIF
 !
@@ -1456 +1460 @@
 !--             Note, each topography is flagged as building. 
                 IF ( zu(k) - ocean_offset <= buildings_f%var_2d(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), 2 ) !indicates building
+                   topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
+                   topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 ) !indicates building
                 ENDIF
              ENDDO
@@ -1712 +1716 @@
     
     USE kinds
+    
+    USE netcdf_data_input_mod,                                                 &
+        ONLY:  buildings_f, terrain_height_f 
 
     USE pegrid
@@ -1749 +1756 @@
     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_local         !< index for topography top at cell-center
     INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  topo !< input array for 3D topography and dummy array for setting "outer"-flags
-
-
+!
+!-- Check for correct setting of the namelist parameter topography. If 
+!-- topography information is read from file but topography = 'flat', 
+!-- initialization does not properly.
+    IF ( ( buildings_f%from_file  .OR.  terrain_height_f%from_file )  .AND.    &
+           TRIM( topography ) /= 'read_from_file' )  THEN
+       message_string =  'If topography information is provided (via ' //      &
+                         'Netcdf or ASCII input) topography = '        //      &
+                         '"read_from_file" is required.'
+       CALL message( 'init_grid', 'PA0437', 1, 2, 0, 6, 0 )      
+    ENDIF
 !
 !-- Set outer and inner index arrays for non-flat topography.
@@ -2263 +2279 @@
     USE control_parameters,                                                    &
         ONLY:  bc_lr_cyc, bc_ns_cyc, constant_flux_layer, land_surface,        &
-               scalar_advec, use_surface_fluxes, use_top_fluxes, urban_surface
+               scalar_advec, topography, use_surface_fluxes, use_top_fluxes,   &
+               urban_surface
 
     USE indices,                                                               &
@@ -2526 +2543 @@
 !
 !-- Finally, set identification flags indicating natural terrain or buildings.
-!-- Natural terrain grid points.
+!-- Natural terrain grid points. Information on the type of the surface is 
+!-- stored in bit 1 of 3D Integer array topo. However, this bit is only set 
+!-- when topography is read from file. In order to run the land-surface model
+!-- also without topography information, set bit 1 explicitely in this case.
     IF ( land_surface )  THEN
-       DO i = nxl, nxr
-          DO j = nys, nyn
-             DO k = nzb, nzt+1
-!
-!--             Natural terrain grid point
-                IF ( BTEST( topo(k,j,i), 1 ) )                                 &
-                   wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 5 )
+!
+!--    If no topography is initialized, the land-surface is at k = nzb.
+       IF ( TRIM( topography ) /= 'read_from_file' )  THEN
+          wall_flags_0(nzb,:,:) = IBSET( wall_flags_0(nzb,:,:), 5 )
+       ELSE
+          DO i = nxl, nxr
+             DO j = nys, nyn
+                DO k = nzb, nzt+1
+!         
+!--                Natural terrain grid point
+                   IF ( BTEST( topo(k,j,i), 1 ) )                              &
+                      wall_flags_0(k,j,i) = IBSET( wall_flags_0(k,j,i), 5 )
+                ENDDO
              ENDDO
           ENDDO
-       ENDDO
+       ENDIF
     ENDIF
 !