Changeset 4102 for palm/trunk/SOURCE/surface_mod.f90
- Timestamp:
- Jul 17, 2019 4:00:03 PM (5 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/surface_mod.f90 (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/surface_mod.f90
r3943 r4102 26 26 ! ----------------- 27 27 ! $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 28 32 ! Removed qsws_eb as it is no longer needed. 29 33 ! … … 306 310 !-- are applied 307 311 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 314 320 315 321 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: start_index !< start index within surface data type for given (j,i) … … 600 606 601 607 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 602 609 603 610 TYPE (surf_type), DIMENSION(0:2), TARGET :: surf_def_h !< horizontal default surfaces (Up, Down, and Top) … … 668 675 ! 669 676 !-- 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, & 672 679 vertical_surfaces_exist, surf_bulk_cloud_model, surf_microphysics_morrison, & 673 680 surf_microphysics_seifert … … 687 694 ! Description: 688 695 ! ------------ 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. 690 698 !------------------------------------------------------------------------------! 691 699 SUBROUTINE init_bc … … 696 704 INTEGER(iwp) :: j !< loop index along y-direction 697 705 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 721 809 ENDDO 722 810 ENDDO 723 811 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 767 831 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 776 870 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 778 876 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) - 1781 start_index_h(0) = bc_h(0)%end_index(j,i) + 1782 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) - 1785 start_index_h(1) = bc_h(1)%end_index(j,i) + 1786 877 ENDDO 787 878 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
