Ignore:
Timestamp:
Jul 17, 2019 4:00:03 PM (2 years ago)
Author:
suehring
Message:

Bugfix, set Neumann boundary conditions for the subgrid TKE at vertical walls instead of implicit Dirichlet conditions that always act as a sink term for the subgrid TKE. Therefore, add new data structure for vertical surfaces and revise the setting of the boundary grid point index space. Moreover, slightly revise setting of boundary conditions at upward- and downward facing surfaces. Finally, setting of boundary conditions for subgrid TKE and dissipation (in RANS mode) is now modularized. Update test case results.

File:
1 edited

Legend:

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

    r3943 r4102  
    2626! -----------------
    2727! $Id$
     28! - Revise initialization of the boundary data structure
     29! - Add new data structure to set boundary conditions at vertical walls
     30!
     31! 3943 2019-05-02 09:50:41Z maronga
    2832! Removed qsws_eb as it is no longer needed.
    2933!
     
    306310!-- are applied
    307311    TYPE bc_type
    308 
    309        INTEGER(iwp) :: ns                                  !< number of surface elements on the PE
    310 
    311        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i       !< x-index linking to the PALM 3D-grid 
    312        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j       !< y-index linking to the PALM 3D-grid   
    313        INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k       !< z-index linking to the PALM 3D-grid   
     312       INTEGER(iwp) :: ioff !< offset value in x-direction, used to determine index of surface element
     313       INTEGER(iwp) :: joff !< offset value in y-direction, used to determine index of surface element
     314       INTEGER(iwp) :: koff !< offset value in z-direction, used to determine index of surface element
     315       INTEGER(iwp) :: ns   !< number of surface elements on the PE
     316
     317       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  i !< x-index linking to the PALM 3D-grid 
     318       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  j !< y-index linking to the PALM 3D-grid   
     319       INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  k !< z-index linking to the PALM 3D-grid   
    314320
    315321       INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: start_index !< start index within surface data type for given (j,i)
     
    600606
    601607    TYPE (bc_type), DIMENSION(0:1)           ::  bc_h        !< boundary condition data type, horizontal upward- and downward facing surfaces
     608    TYPE (bc_type), DIMENSION(0:3)           ::  bc_v        !< boundary condition data type, vertical surfaces
    602609
    603610    TYPE (surf_type), DIMENSION(0:2), TARGET ::  surf_def_h  !< horizontal default surfaces (Up, Down, and Top)
     
    668675!
    669676!-- Public variables
    670     PUBLIC bc_h, ind_pav_green, ind_veg_wall, ind_wat_win, ns_h_on_file, ns_v_on_file, surf_def_h, &
    671            surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v, surf_type,                  &
     677    PUBLIC bc_h, bc_v, ind_pav_green, ind_veg_wall, ind_wat_win, ns_h_on_file, ns_v_on_file,      &
     678           surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v, surf_type,      &
    672679           vertical_surfaces_exist, surf_bulk_cloud_model, surf_microphysics_morrison,             &
    673680           surf_microphysics_seifert
     
    687694! Description:
    688695! ------------
    689 !> Initialize data type for setting boundary conditions at horizontal surfaces.
     696!> Initialize data type for setting boundary conditions at horizontal and 
     697!> vertical surfaces.
    690698!------------------------------------------------------------------------------!
    691699    SUBROUTINE init_bc
     
    696704       INTEGER(iwp) ::  j         !< loop index along y-direction
    697705       INTEGER(iwp) ::  k         !< loop index along y-direction
    698 
    699        INTEGER(iwp), DIMENSION(0:1) ::  num_h         !< number of surfaces
    700        INTEGER(iwp), DIMENSION(0:1) ::  num_h_kji     !< number of surfaces for (j,i)-grid point
    701        INTEGER(iwp), DIMENSION(0:1) ::  start_index_h !< local start index
    702 
    703 !
    704 !--    First of all, count the number of upward- and downward-facing surfaces
    705        num_h = 0
    706        DO  i = nxlg, nxrg
    707           DO  j = nysg, nyng
    708              DO  k = nzb+1, nzt
    709 !
    710 !--             Check if current gridpoint belongs to the atmosphere
    711                 IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
    712 !
    713 !--                Upward-facing
    714                    IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) )              &
    715                       num_h(0) = num_h(0) + 1
    716 !
    717 !--                Downward-facing
    718                    IF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) )              &
    719                       num_h(1) = num_h(1) + 1
    720                 ENDIF
     706       INTEGER(iwp) ::  l         !< running index for differently aligned surfaces
     707
     708       INTEGER(iwp), DIMENSION(0:1) ::  num_h         !< number of horizontal surfaces on subdomain
     709       INTEGER(iwp), DIMENSION(0:1) ::  num_h_kji     !< number of horizontal surfaces at (j,i)-grid point
     710       INTEGER(iwp), DIMENSION(0:1) ::  start_index_h !< local start index of horizontal surface elements
     711       
     712       INTEGER(iwp), DIMENSION(0:3) ::  num_v         !< number of vertical surfaces on subdomain
     713       INTEGER(iwp), DIMENSION(0:3) ::  num_v_kji     !< number of vertical surfaces at (j,i)-grid point
     714       INTEGER(iwp), DIMENSION(0:3) ::  start_index_v !< local start index of vertical surface elements
     715!
     716!--    Set offset indices, i.e. index difference between surface element and
     717!--    surface-bounded grid point.
     718!--    Horizontal surfaces - no horizontal offsets
     719       bc_h(:)%ioff = 0
     720       bc_h(:)%joff = 0
     721!
     722!--    Horizontal surfaces, upward facing (0) and downward facing (1)
     723       bc_h(0)%koff   = -1
     724       bc_h(1)%koff   = 1
     725!
     726!--    Vertical surfaces - no vertical offset
     727       bc_v(0:3)%koff = 0
     728!
     729!--    North- and southward facing - no offset in x
     730       bc_v(0:1)%ioff = 0
     731!
     732!--    Northward facing offset in y
     733       bc_v(0)%joff = -1
     734!
     735!--    Southward facing offset in y
     736       bc_v(1)%joff = 1
     737!
     738!--    East- and westward facing - no offset in y
     739       bc_v(2:3)%joff = 0
     740!
     741!--    Eastward facing offset in x
     742       bc_v(2)%ioff = -1
     743!
     744!--    Westward facing offset in y
     745       bc_v(3)%ioff = 1
     746!
     747!--    Initialize data structure for horizontal surfaces, i.e. count the number
     748!--    of surface elements, allocate and initialize the respective index arrays,
     749!--    and set the respective start and end indices at each (j,i)-location.
     750       DO  l = 0, 1
     751!
     752!--       Count the number of upward- and downward-facing surfaces on subdomain
     753          num_h(l) = 0
     754          DO  i = nxlg, nxrg
     755             DO  j = nysg, nyng
     756                DO  k = nzb+1, nzt
     757!         
     758!--                Check if current gridpoint belongs to the atmosphere
     759                   IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
     760                      IF ( .NOT. BTEST( wall_flags_0(k+bc_h(l)%koff,           &
     761                                                     j+bc_h(l)%joff,           &
     762                                                     i+bc_h(l)%ioff), 0 ) )    &
     763                         num_h(l) = num_h(l) + 1
     764                   ENDIF
     765                ENDDO
     766             ENDDO
     767          ENDDO
     768!         
     769!--       Save the number of horizontal surface elements
     770          bc_h(l)%ns = num_h(l)
     771!
     772!--       ALLOCATE arrays for horizontal surfaces
     773          ALLOCATE( bc_h(l)%i(1:bc_h(l)%ns) )
     774          ALLOCATE( bc_h(l)%j(1:bc_h(l)%ns) )
     775          ALLOCATE( bc_h(l)%k(1:bc_h(l)%ns) )
     776          ALLOCATE( bc_h(l)%start_index(nysg:nyng,nxlg:nxrg) )
     777          ALLOCATE( bc_h(l)%end_index(nysg:nyng,nxlg:nxrg)   )
     778          bc_h(l)%start_index = 1
     779          bc_h(l)%end_index   = 0
     780         
     781          num_h(l)         = 1
     782          start_index_h(l) = 1
     783          DO  i = nxlg, nxrg
     784             DO  j = nysg, nyng
     785         
     786                num_h_kji(l) = 0
     787                DO  k = nzb+1, nzt
     788!         
     789!--                Check if current gridpoint belongs to the atmosphere
     790                   IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
     791!         
     792!--                   Upward-facing
     793                      IF ( .NOT. BTEST( wall_flags_0(k+bc_h(l)%koff,           &
     794                                                     j+bc_h(l)%joff,           &
     795                                                     i+bc_h(l)%ioff), 0 )      &
     796                         )  THEN
     797                         bc_h(l)%i(num_h(l)) = i
     798                         bc_h(l)%j(num_h(l)) = j
     799                         bc_h(l)%k(num_h(l)) = k
     800                         num_h_kji(l)        = num_h_kji(l) + 1
     801                         num_h(l)            = num_h(l) + 1
     802                      ENDIF
     803                   ENDIF
     804                ENDDO
     805                bc_h(l)%start_index(j,i) = start_index_h(l)
     806                bc_h(l)%end_index(j,i)   = bc_h(l)%start_index(j,i) +          &
     807                                           num_h_kji(l) - 1
     808                start_index_h(l)         = bc_h(l)%end_index(j,i) + 1
    721809             ENDDO
    722810          ENDDO
    723811       ENDDO
    724 !
    725 !--    Save the number of surface elements
    726        bc_h(0)%ns = num_h(0)
    727        bc_h(1)%ns = num_h(1)
    728 !
    729 !--    ALLOCATE data type variables
    730 !--    Upward facing
    731        ALLOCATE( bc_h(0)%i(1:bc_h(0)%ns) )
    732        ALLOCATE( bc_h(0)%j(1:bc_h(0)%ns) )
    733        ALLOCATE( bc_h(0)%k(1:bc_h(0)%ns) )
    734        ALLOCATE( bc_h(0)%start_index(nysg:nyng,nxlg:nxrg) )
    735        ALLOCATE( bc_h(0)%end_index(nysg:nyng,nxlg:nxrg)   )
    736        bc_h(0)%start_index = 1
    737        bc_h(0)%end_index   = 0
    738 !
    739 !--    Downward facing
    740        ALLOCATE( bc_h(1)%i(1:bc_h(1)%ns) )
    741        ALLOCATE( bc_h(1)%j(1:bc_h(1)%ns) )
    742        ALLOCATE( bc_h(1)%k(1:bc_h(1)%ns) )
    743        ALLOCATE( bc_h(1)%start_index(nysg:nyng,nxlg:nxrg) )
    744        ALLOCATE( bc_h(1)%end_index(nysg:nyng,nxlg:nxrg)   )
    745        bc_h(1)%start_index = 1
    746        bc_h(1)%end_index   = 0
    747 !
    748 !--    Store the respective indices on data type
    749        num_h(0:1)         = 1
    750        start_index_h(0:1) = 1
    751        DO  i = nxlg, nxrg
    752           DO  j = nysg, nyng
    753 
    754              num_h_kji(0:1) = 0
    755              DO  k = nzb+1, nzt
    756 !
    757 !--             Check if current gridpoint belongs to the atmosphere
    758                 IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
    759 !
    760 !--                Upward-facing
    761                    IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) )  THEN
    762                       bc_h(0)%i(num_h(0)) = i
    763                       bc_h(0)%j(num_h(0)) = j
    764                       bc_h(0)%k(num_h(0)) = k
    765                       num_h_kji(0)        = num_h_kji(0) + 1
    766                       num_h(0)            = num_h(0) + 1
     812
     813!
     814!--    Initialize data structure for vertical surfaces, i.e. count the number
     815!--    of surface elements, allocate and initialize the respective index arrays,
     816!--    and set the respective start and end indices at each (j,i)-location.
     817       DO  l = 0, 3
     818!
     819!--       Count the number of upward- and downward-facing surfaces on subdomain
     820          num_v(l) = 0
     821          DO  i = nxlg, nxrg
     822             DO  j = nysg, nyng
     823                DO  k = nzb+1, nzt
     824!         
     825!--                Check if current gridpoint belongs to the atmosphere
     826                   IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
     827                      IF ( .NOT. BTEST( wall_flags_0(k+bc_v(l)%koff,           &
     828                                                     j+bc_v(l)%joff,           &
     829                                                     i+bc_v(l)%ioff), 0 ) )    &
     830                         num_v(l) = num_v(l) + 1
    767831                   ENDIF
    768 !
    769 !--                Downward-facing
    770                    IF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) )  THEN
    771                       bc_h(1)%i(num_h(1)) = i
    772                       bc_h(1)%j(num_h(1)) = j
    773                       bc_h(1)%k(num_h(1)) = k
    774                       num_h_kji(1)        = num_h_kji(1) + 1
    775                       num_h(1)            = num_h(1) + 1
     832                ENDDO
     833             ENDDO
     834          ENDDO
     835!         
     836!--       Save the number of horizontal surface elements
     837          bc_v(l)%ns = num_v(l)
     838!
     839!--       ALLOCATE arrays for horizontal surfaces
     840          ALLOCATE( bc_v(l)%i(1:bc_v(l)%ns) )
     841          ALLOCATE( bc_v(l)%j(1:bc_v(l)%ns) )
     842          ALLOCATE( bc_v(l)%k(1:bc_v(l)%ns) )
     843          ALLOCATE( bc_v(l)%start_index(nysg:nyng,nxlg:nxrg) )
     844          ALLOCATE( bc_v(l)%end_index(nysg:nyng,nxlg:nxrg)   )
     845          bc_v(l)%start_index = 1
     846          bc_v(l)%end_index   = 0
     847         
     848          num_v(l)         = 1
     849          start_index_v(l) = 1
     850          DO  i = nxlg, nxrg
     851             DO  j = nysg, nyng
     852         
     853                num_v_kji(l) = 0
     854                DO  k = nzb+1, nzt
     855!         
     856!--                Check if current gridpoint belongs to the atmosphere
     857                   IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
     858!         
     859!--                   Upward-facing
     860                      IF ( .NOT. BTEST( wall_flags_0(k+bc_v(l)%koff,           &
     861                                                     j+bc_v(l)%joff,           &
     862                                                     i+bc_v(l)%ioff), 0 )      &
     863                         )  THEN
     864                         bc_v(l)%i(num_v(l)) = i
     865                         bc_v(l)%j(num_v(l)) = j
     866                         bc_v(l)%k(num_v(l)) = k
     867                         num_v_kji(l)        = num_v_kji(l) + 1
     868                         num_v(l)            = num_v(l) + 1
     869                      ENDIF
    776870                   ENDIF
    777                 ENDIF
     871                ENDDO
     872                bc_v(l)%start_index(j,i) = start_index_v(l)
     873                bc_v(l)%end_index(j,i)   = bc_v(l)%start_index(j,i) +          &
     874                                           num_v_kji(l) - 1
     875                start_index_v(l)         = bc_v(l)%end_index(j,i) + 1
    778876             ENDDO
    779              bc_h(0)%start_index(j,i) = start_index_h(0)
    780              bc_h(0)%end_index(j,i)   = bc_h(0)%start_index(j,i) + num_h_kji(0) - 1
    781              start_index_h(0)         = bc_h(0)%end_index(j,i) + 1
    782 
    783              bc_h(1)%start_index(j,i) = start_index_h(1)
    784              bc_h(1)%end_index(j,i)   = bc_h(1)%start_index(j,i) + num_h_kji(1) - 1
    785              start_index_h(1)         = bc_h(1)%end_index(j,i) + 1
    786877          ENDDO
    787878       ENDDO
Note: See TracChangeset for help on using the changeset viewer.