Index: palm/trunk/SOURCE/surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/surface_mod.f90	(revision 3943)
+++ palm/trunk/SOURCE/surface_mod.f90	(revision 4102)
@@ -26,4 +26,8 @@
 ! -----------------
 ! $Id$
+! - Revise initialization of the boundary data structure
+! - Add new data structure to set boundary conditions at vertical walls
+! 
+! 3943 2019-05-02 09:50:41Z maronga
 ! Removed qsws_eb as it is no longer needed.
 ! 
@@ -306,10 +310,12 @@
 !-- are applied 
     TYPE bc_type
-
-       INTEGER(iwp) :: ns                                  !< number of surface elements on the PE
-
-       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i       !< x-index linking to the PALM 3D-grid  
-       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j       !< y-index linking to the PALM 3D-grid    
-       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k       !< z-index linking to the PALM 3D-grid   
+       INTEGER(iwp) :: ioff !< offset value in x-direction, used to determine index of surface element
+       INTEGER(iwp) :: joff !< offset value in y-direction, used to determine index of surface element
+       INTEGER(iwp) :: koff !< offset value in z-direction, used to determine index of surface element
+       INTEGER(iwp) :: ns   !< number of surface elements on the PE
+
+       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i !< x-index linking to the PALM 3D-grid  
+       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j !< y-index linking to the PALM 3D-grid    
+       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k !< z-index linking to the PALM 3D-grid   
 
        INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: start_index !< start index within surface data type for given (j,i) 
@@ -600,4 +606,5 @@
 
     TYPE (bc_type), DIMENSION(0:1)           ::  bc_h        !< boundary condition data type, horizontal upward- and downward facing surfaces
+    TYPE (bc_type), DIMENSION(0:3)           ::  bc_v        !< boundary condition data type, vertical surfaces
 
     TYPE (surf_type), DIMENSION(0:2), TARGET ::  surf_def_h  !< horizontal default surfaces (Up, Down, and Top)
@@ -668,6 +675,6 @@
 !
 !-- Public variables
-    PUBLIC bc_h, ind_pav_green, ind_veg_wall, ind_wat_win, ns_h_on_file, ns_v_on_file, surf_def_h, &
-           surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v, surf_type,                  &
+    PUBLIC bc_h, bc_v, ind_pav_green, ind_veg_wall, ind_wat_win, ns_h_on_file, ns_v_on_file,       &
+           surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v, surf_type,      &
            vertical_surfaces_exist, surf_bulk_cloud_model, surf_microphysics_morrison,             &
            surf_microphysics_seifert
@@ -687,5 +694,6 @@
 ! Description:
 ! ------------
-!> Initialize data type for setting boundary conditions at horizontal surfaces. 
+!> Initialize data type for setting boundary conditions at horizontal and  
+!> vertical surfaces. 
 !------------------------------------------------------------------------------!
     SUBROUTINE init_bc
@@ -696,92 +704,175 @@
        INTEGER(iwp) ::  j         !< loop index along y-direction
        INTEGER(iwp) ::  k         !< loop index along y-direction
-
-       INTEGER(iwp), DIMENSION(0:1) ::  num_h         !< number of surfaces
-       INTEGER(iwp), DIMENSION(0:1) ::  num_h_kji     !< number of surfaces for (j,i)-grid point
-       INTEGER(iwp), DIMENSION(0:1) ::  start_index_h !< local start index
-
-!
-!--    First of all, count the number of upward- and downward-facing surfaces
-       num_h = 0
-       DO  i = nxlg, nxrg
-          DO  j = nysg, nyng
-             DO  k = nzb+1, nzt
-!
-!--             Check if current gridpoint belongs to the atmosphere
-                IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
-!
-!--                Upward-facing
-                   IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) )              &
-                      num_h(0) = num_h(0) + 1
-!
-!--                Downward-facing
-                   IF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) )              &
-                      num_h(1) = num_h(1) + 1
-                ENDIF
+       INTEGER(iwp) ::  l         !< running index for differently aligned surfaces
+
+       INTEGER(iwp), DIMENSION(0:1) ::  num_h         !< number of horizontal surfaces on subdomain
+       INTEGER(iwp), DIMENSION(0:1) ::  num_h_kji     !< number of horizontal surfaces at (j,i)-grid point
+       INTEGER(iwp), DIMENSION(0:1) ::  start_index_h !< local start index of horizontal surface elements
+       
+       INTEGER(iwp), DIMENSION(0:3) ::  num_v         !< number of vertical surfaces on subdomain
+       INTEGER(iwp), DIMENSION(0:3) ::  num_v_kji     !< number of vertical surfaces at (j,i)-grid point
+       INTEGER(iwp), DIMENSION(0:3) ::  start_index_v !< local start index of vertical surface elements
+!
+!--    Set offset indices, i.e. index difference between surface element and 
+!--    surface-bounded grid point.
+!--    Horizontal surfaces - no horizontal offsets
+       bc_h(:)%ioff = 0
+       bc_h(:)%joff = 0
+!
+!--    Horizontal surfaces, upward facing (0) and downward facing (1)
+       bc_h(0)%koff   = -1
+       bc_h(1)%koff   = 1
+!
+!--    Vertical surfaces - no vertical offset
+       bc_v(0:3)%koff = 0
+!
+!--    North- and southward facing - no offset in x
+       bc_v(0:1)%ioff = 0
+!
+!--    Northward facing offset in y
+       bc_v(0)%joff = -1
+!
+!--    Southward facing offset in y
+       bc_v(1)%joff = 1
+!
+!--    East- and westward facing - no offset in y
+       bc_v(2:3)%joff = 0
+!
+!--    Eastward facing offset in x
+       bc_v(2)%ioff = -1
+!
+!--    Westward facing offset in y
+       bc_v(3)%ioff = 1
+!
+!--    Initialize data structure for horizontal surfaces, i.e. count the number
+!--    of surface elements, allocate and initialize the respective index arrays, 
+!--    and set the respective start and end indices at each (j,i)-location.
+       DO  l = 0, 1
+!
+!--       Count the number of upward- and downward-facing surfaces on subdomain
+          num_h(l) = 0
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
+                DO  k = nzb+1, nzt
+!         
+!--                Check if current gridpoint belongs to the atmosphere
+                   IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
+                      IF ( .NOT. BTEST( wall_flags_0(k+bc_h(l)%koff,           &
+                                                     j+bc_h(l)%joff,           &
+                                                     i+bc_h(l)%ioff), 0 ) )    &
+                         num_h(l) = num_h(l) + 1
+                   ENDIF
+                ENDDO
+             ENDDO
+          ENDDO
+!         
+!--       Save the number of horizontal surface elements
+          bc_h(l)%ns = num_h(l)
+!
+!--       ALLOCATE arrays for horizontal surfaces
+          ALLOCATE( bc_h(l)%i(1:bc_h(l)%ns) )
+          ALLOCATE( bc_h(l)%j(1:bc_h(l)%ns) )
+          ALLOCATE( bc_h(l)%k(1:bc_h(l)%ns) )
+          ALLOCATE( bc_h(l)%start_index(nysg:nyng,nxlg:nxrg) )
+          ALLOCATE( bc_h(l)%end_index(nysg:nyng,nxlg:nxrg)   )
+          bc_h(l)%start_index = 1
+          bc_h(l)%end_index   = 0
+          
+          num_h(l)         = 1
+          start_index_h(l) = 1
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
+          
+                num_h_kji(l) = 0
+                DO  k = nzb+1, nzt
+!         
+!--                Check if current gridpoint belongs to the atmosphere
+                   IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
+!         
+!--                   Upward-facing
+                      IF ( .NOT. BTEST( wall_flags_0(k+bc_h(l)%koff,           &
+                                                     j+bc_h(l)%joff,           &
+                                                     i+bc_h(l)%ioff), 0 )      &
+                         )  THEN
+                         bc_h(l)%i(num_h(l)) = i
+                         bc_h(l)%j(num_h(l)) = j
+                         bc_h(l)%k(num_h(l)) = k
+                         num_h_kji(l)        = num_h_kji(l) + 1
+                         num_h(l)            = num_h(l) + 1
+                      ENDIF
+                   ENDIF
+                ENDDO
+                bc_h(l)%start_index(j,i) = start_index_h(l)
+                bc_h(l)%end_index(j,i)   = bc_h(l)%start_index(j,i) +          &
+                                           num_h_kji(l) - 1
+                start_index_h(l)         = bc_h(l)%end_index(j,i) + 1
              ENDDO
           ENDDO
        ENDDO
-!
-!--    Save the number of surface elements
-       bc_h(0)%ns = num_h(0)
-       bc_h(1)%ns = num_h(1)
-!
-!--    ALLOCATE data type variables
-!--    Upward facing
-       ALLOCATE( bc_h(0)%i(1:bc_h(0)%ns) )
-       ALLOCATE( bc_h(0)%j(1:bc_h(0)%ns) )
-       ALLOCATE( bc_h(0)%k(1:bc_h(0)%ns) )
-       ALLOCATE( bc_h(0)%start_index(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE( bc_h(0)%end_index(nysg:nyng,nxlg:nxrg)   )
-       bc_h(0)%start_index = 1
-       bc_h(0)%end_index   = 0
-!
-!--    Downward facing
-       ALLOCATE( bc_h(1)%i(1:bc_h(1)%ns) )
-       ALLOCATE( bc_h(1)%j(1:bc_h(1)%ns) )
-       ALLOCATE( bc_h(1)%k(1:bc_h(1)%ns) )
-       ALLOCATE( bc_h(1)%start_index(nysg:nyng,nxlg:nxrg) )
-       ALLOCATE( bc_h(1)%end_index(nysg:nyng,nxlg:nxrg)   )
-       bc_h(1)%start_index = 1
-       bc_h(1)%end_index   = 0
-!
-!--    Store the respective indices on data type
-       num_h(0:1)         = 1
-       start_index_h(0:1) = 1
-       DO  i = nxlg, nxrg
-          DO  j = nysg, nyng
-
-             num_h_kji(0:1) = 0
-             DO  k = nzb+1, nzt
-!
-!--             Check if current gridpoint belongs to the atmosphere
-                IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
-!
-!--                Upward-facing
-                   IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) )  THEN
-                      bc_h(0)%i(num_h(0)) = i
-                      bc_h(0)%j(num_h(0)) = j
-                      bc_h(0)%k(num_h(0)) = k
-                      num_h_kji(0)        = num_h_kji(0) + 1
-                      num_h(0)            = num_h(0) + 1
+
+!
+!--    Initialize data structure for vertical surfaces, i.e. count the number
+!--    of surface elements, allocate and initialize the respective index arrays, 
+!--    and set the respective start and end indices at each (j,i)-location.
+       DO  l = 0, 3
+!
+!--       Count the number of upward- and downward-facing surfaces on subdomain
+          num_v(l) = 0
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
+                DO  k = nzb+1, nzt
+!         
+!--                Check if current gridpoint belongs to the atmosphere
+                   IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
+                      IF ( .NOT. BTEST( wall_flags_0(k+bc_v(l)%koff,           &
+                                                     j+bc_v(l)%joff,           &
+                                                     i+bc_v(l)%ioff), 0 ) )    &
+                         num_v(l) = num_v(l) + 1
                    ENDIF
-!
-!--                Downward-facing
-                   IF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) )  THEN
-                      bc_h(1)%i(num_h(1)) = i
-                      bc_h(1)%j(num_h(1)) = j
-                      bc_h(1)%k(num_h(1)) = k
-                      num_h_kji(1)        = num_h_kji(1) + 1
-                      num_h(1)            = num_h(1) + 1
+                ENDDO
+             ENDDO
+          ENDDO
+!         
+!--       Save the number of horizontal surface elements
+          bc_v(l)%ns = num_v(l)
+!
+!--       ALLOCATE arrays for horizontal surfaces
+          ALLOCATE( bc_v(l)%i(1:bc_v(l)%ns) )
+          ALLOCATE( bc_v(l)%j(1:bc_v(l)%ns) )
+          ALLOCATE( bc_v(l)%k(1:bc_v(l)%ns) )
+          ALLOCATE( bc_v(l)%start_index(nysg:nyng,nxlg:nxrg) )
+          ALLOCATE( bc_v(l)%end_index(nysg:nyng,nxlg:nxrg)   )
+          bc_v(l)%start_index = 1
+          bc_v(l)%end_index   = 0
+          
+          num_v(l)         = 1
+          start_index_v(l) = 1
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
+          
+                num_v_kji(l) = 0
+                DO  k = nzb+1, nzt
+!         
+!--                Check if current gridpoint belongs to the atmosphere
+                   IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
+!         
+!--                   Upward-facing
+                      IF ( .NOT. BTEST( wall_flags_0(k+bc_v(l)%koff,           &
+                                                     j+bc_v(l)%joff,           &
+                                                     i+bc_v(l)%ioff), 0 )      &
+                         )  THEN
+                         bc_v(l)%i(num_v(l)) = i
+                         bc_v(l)%j(num_v(l)) = j
+                         bc_v(l)%k(num_v(l)) = k
+                         num_v_kji(l)        = num_v_kji(l) + 1
+                         num_v(l)            = num_v(l) + 1
+                      ENDIF
                    ENDIF
-                ENDIF
+                ENDDO
+                bc_v(l)%start_index(j,i) = start_index_v(l)
+                bc_v(l)%end_index(j,i)   = bc_v(l)%start_index(j,i) +          &
+                                           num_v_kji(l) - 1
+                start_index_v(l)         = bc_v(l)%end_index(j,i) + 1
              ENDDO
-             bc_h(0)%start_index(j,i) = start_index_h(0)
-             bc_h(0)%end_index(j,i)   = bc_h(0)%start_index(j,i) + num_h_kji(0) - 1
-             start_index_h(0)         = bc_h(0)%end_index(j,i) + 1
-
-             bc_h(1)%start_index(j,i) = start_index_h(1)
-             bc_h(1)%end_index(j,i)   = bc_h(1)%start_index(j,i) + num_h_kji(1) - 1
-             start_index_h(1)         = bc_h(1)%end_index(j,i) + 1
           ENDDO
        ENDDO
