Changeset 4414 for palm/trunk/SOURCE/init_grid.f90
- Timestamp:
- Feb 19, 2020 8:16:04 PM (4 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/init_grid.f90 (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/init_grid.f90
r4386 r4414 25 25 ! ----------------- 26 26 ! $Id$ 27 ! - Remove deprecated topography arrays nzb_s_inner, nzb_u_inner, etc. 28 ! - Move initialization of boundary conditions and multigrid into an extra 29 ! module interface. 30 ! 31 ! 4386 2020-01-27 15:07:30Z Giersch 27 32 ! Allocation statements, comments, naming of variables revised and _wp added to 28 33 ! real type values … … 123 128 !------------------------------------------------------------------------------! 124 129 SUBROUTINE init_grid 125 126 USE advec_ws, &127 ONLY: ws_init_flags_momentum, ws_init_flags_scalar128 130 129 131 USE arrays_3d, & 130 132 ONLY: dd2zu, ddzu, ddzu_pres, ddzw, dzu, dzw, x, xu, y, yv, zu, zw 131 133 132 134 USE control_parameters, & 133 ONLY: bc_lr_cyc, bc_ns_cyc, & 134 bc_dirichlet_l, & 135 bc_dirichlet_n, & 136 bc_dirichlet_r, & 137 bc_dirichlet_s, & 138 bc_radiation_l, & 139 bc_radiation_n, & 140 bc_radiation_r, & 141 bc_radiation_s, & 142 constant_flux_layer, dz, dz_max, dz_stretch_factor, & 135 ONLY: constant_flux_layer, dz, dz_max, dz_stretch_factor, & 143 136 dz_stretch_factor_array, dz_stretch_level, dz_stretch_level_end,& 144 137 dz_stretch_level_end_index, dz_stretch_level_start_index, & 145 138 dz_stretch_level_start, ibc_uv_b, message_string, & 146 momentum_advec, number_stretch_level_end, & 147 number_stretch_level_start, ocean_mode, psolver, scalar_advec, & 148 symmetry_flag, topography, use_surface_fluxes 149 139 number_stretch_level_end, & 140 number_stretch_level_start, & 141 ocean_mode, & 142 psolver, & 143 symmetry_flag, & 144 topography, & 145 use_surface_fluxes 146 150 147 USE grid_variables, & 151 148 ONLY: ddx, ddx2, ddy, ddy2, dx, dx2, dy, dy2, zu_s_inner, zw_w_inner 152 149 153 150 USE indices, & 154 ONLY: advc_flags_m, & 155 advc_flags_s, & 156 nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz, & 157 nzb, nzb_diff, nzb_diff_s_inner, nzb_diff_s_outer, & 158 nzb_max, nzb_s_inner, nzb_s_outer, nzb_u_inner, & 159 nzb_u_outer, nzb_v_inner, nzb_v_outer, nzb_w_inner, & 160 nzb_w_outer, nzt, topo_top_ind, topo_min_level 161 151 ONLY: nbgp, & 152 nx, & 153 nxl, & 154 nxlg, & 155 nxr, & 156 nxrg, & 157 ny, & 158 nyn, & 159 nyng, & 160 nys, & 161 nysg, & 162 nz, & 163 nzb, & 164 nzb_diff, & 165 nzb_max, & 166 nzt, & 167 topo_top_ind, & 168 topo_min_level 169 162 170 USE kinds 163 171 164 172 USE pegrid 165 166 USE poismg_noopt_mod167 168 USE surface_mod, &169 ONLY: init_bc170 173 171 174 USE vertical_nesting_mod, & … … 182 185 INTEGER(iwp) :: nzb_local_max !< vertical grid index of maximum topography height 183 186 INTEGER(iwp) :: nzb_local_min !< vertical grid index of minimum topography height 184 185 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nzb_local !< index for topography top at cell-center186 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nzb_tmp !< dummy to calculate topography indices on u- and v-grid187 187 188 188 INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: topo !< input array for 3D topography and dummy array for setting "outer"-flags … … 644 644 ! 645 645 !-- Set flags to mask topography on the grid. 646 CALL set_topo_flags( topo ) 647 ! 648 !-- Calculate wall flag arrays for the multigrid method. 649 !-- Please note, wall flags are only applied in the non-optimized version. 650 IF ( psolver == 'multigrid_noopt' ) CALL poismg_noopt_init 651 652 ! 653 !-- Init flags for ws-scheme to degrade order of the numerics near walls, i.e. 654 !-- to decrease the numerical stencil appropriately. The order of the scheme 655 !-- is degraded near solid walls as well as near non-cyclic inflow and outflow 656 !-- boundaries. Do this separately for momentum and scalars. 657 IF ( momentum_advec == 'ws-scheme' ) THEN 658 ALLOCATE( advc_flags_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 659 CALL ws_init_flags_momentum 660 ENDIF 661 IF ( scalar_advec == 'ws-scheme' ) THEN 662 ALLOCATE( advc_flags_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 663 advc_flags_s = 0 664 665 CALL ws_init_flags_scalar( bc_dirichlet_l .OR. bc_radiation_l, & 666 bc_dirichlet_n .OR. bc_radiation_n, & 667 bc_dirichlet_r .OR. bc_radiation_r, & 668 bc_dirichlet_s .OR. bc_radiation_s, & 669 advc_flags_s ) 670 ENDIF 646 CALL set_topo_flags( topo ) 671 647 672 648 ! … … 704 680 topo_min_level = MINVAL( topo_top_ind(nys:nyn,nxl:nxr,0) ) 705 681 #endif 706 ! 707 !-- Initialize boundary conditions via surface type 708 CALL init_bc 709 710 ! 711 !-- Allocate and set topography height arrays required for data output 712 IF ( TRIM( topography ) /= 'flat' ) THEN 713 ! 714 !-- Allocate and set the arrays containing the topography height 715 IF ( nxr == nx .AND. nyn /= ny ) THEN 716 ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn) ) 717 ALLOCATE( zw_w_inner(nxl:nxr+1,nys:nyn) ) 718 ELSEIF ( nxr /= nx .AND. nyn == ny ) THEN 719 ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn+1) ) 720 ALLOCATE( zw_w_inner(nxl:nxr,nys:nyn+1) ) 721 ELSEIF ( nxr == nx .AND. nyn == ny ) THEN 722 ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn+1) ) 723 ALLOCATE( zw_w_inner(nxl:nxr+1,nys:nyn+1) ) 724 ELSE 725 ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn) ) 726 ALLOCATE( zw_w_inner(nxl:nxr,nys:nyn) ) 727 ENDIF 728 729 zu_s_inner = 0.0_wp 730 zw_w_inner = 0.0_wp 731 ! 732 !-- Determine local topography height on scalar and w-grid. Note, setting 733 !-- lateral boundary values is not necessary, realized via wall_flags_static_0 734 !-- array. Further, please note that loop bounds are different from 735 !-- nxl to nxr and nys to nyn on south and right model boundary, hence, 736 !-- use intrinsic lbound and ubound functions to infer array bounds. 737 DO i = LBOUND(zu_s_inner, 1), UBOUND(zu_s_inner, 1) 738 DO j = LBOUND(zu_s_inner, 2), UBOUND(zu_s_inner, 2) 739 ! 740 !-- Topography height on scalar grid. Therefore, determine index of 741 !-- upward-facing surface element on scalar grid. 742 zu_s_inner(i,j) = zu(topo_top_ind(j,i,0)) 743 ! 744 !-- Topography height on w grid. Therefore, determine index of 745 !-- upward-facing surface element on w grid. 746 zw_w_inner(i,j) = zw(topo_top_ind(j,i,3)) 747 ENDDO 748 ENDDO 749 ENDIF 750 751 ! 752 !-- In the following, calculate 2D index arrays. Note, these will be removed 753 !-- soon. 754 !-- Allocate outer and inner index arrays for topography and set 755 !-- defaults. 756 ALLOCATE( nzb_s_inner(nysg:nyng,nxlg:nxrg) ) 757 ALLOCATE( nzb_s_outer(nysg:nyng,nxlg:nxrg) ) 758 ALLOCATE( nzb_u_inner(nysg:nyng,nxlg:nxrg) ) 759 ALLOCATE( nzb_u_outer(nysg:nyng,nxlg:nxrg) ) 760 ALLOCATE( nzb_v_inner(nysg:nyng,nxlg:nxrg) ) 761 ALLOCATE( nzb_v_outer(nysg:nyng,nxlg:nxrg) ) 762 ALLOCATE( nzb_w_inner(nysg:nyng,nxlg:nxrg) ) 763 ALLOCATE( nzb_w_outer(nysg:nyng,nxlg:nxrg) ) 764 ALLOCATE( nzb_diff_s_inner(nysg:nyng,nxlg:nxrg) ) 765 ALLOCATE( nzb_diff_s_outer(nysg:nyng,nxlg:nxrg) ) 766 ALLOCATE( nzb_local(nysg:nyng,nxlg:nxrg) ) 767 ALLOCATE( nzb_tmp(nysg:nyng,nxlg:nxrg) ) 768 ! 769 !-- Initialize 2D-index arrays. Note, these will be removed soon! 770 nzb_local(nys:nyn,nxl:nxr) = topo_top_ind(nys:nyn,nxl:nxr,0) 771 CALL exchange_horiz_2d_int( nzb_local, nys, nyn, nxl, nxr, nbgp ) 682 772 683 ! 773 684 !-- Check topography for consistency with model domain. Therefore, use … … 792 703 ENDIF 793 704 ENDIF 794 795 nzb_s_inner = nzb; nzb_s_outer = nzb796 nzb_u_inner = nzb; nzb_u_outer = nzb797 nzb_v_inner = nzb; nzb_v_outer = nzb798 nzb_w_inner = nzb; nzb_w_outer = nzb799 800 705 ! 801 706 !-- Define vertical gridpoint from (or to) which on the usual finite difference … … 807 712 ENDIF 808 713 809 nzb_diff_s_inner = nzb_diff; nzb_diff_s_outer = nzb_diff 810 ! 811 !-- Set Neumann conditions for topography. Will be removed soon. 812 IF ( .NOT. bc_ns_cyc ) THEN 813 IF ( nys == 0 ) THEN 814 DO i = 1, nbgp 815 nzb_local(nys-i,:) = nzb_local(nys,:) 816 ENDDO 817 ELSEIF ( nyn == ny ) THEN 818 DO i = 1, nbgp 819 nzb_local(ny+i,:) = nzb_local(ny,:) 820 ENDDO 821 ENDIF 822 ENDIF 823 824 IF ( .NOT. bc_lr_cyc ) THEN 825 IF ( nxl == 0 ) THEN 826 DO i = 1, nbgp 827 nzb_local(:,nxl-i) = nzb_local(:,nxl) 828 ENDDO 829 ELSEIF ( nxr == nx ) THEN 830 DO i = 1, nbgp 831 nzb_local(:,nx+i) = nzb_local(:,nx) 832 ENDDO 833 ENDIF 834 ENDIF 835 ! 836 !-- Initialization of 2D index arrays, will be removed soon! 837 !-- Initialize nzb_s_inner and nzb_w_inner 838 nzb_s_inner = nzb_local 839 nzb_w_inner = nzb_local 840 841 ! 842 !-- Initialize remaining index arrays: 843 !-- first pre-initialize them with nzb_s_inner... 844 nzb_u_inner = nzb_s_inner 845 nzb_u_outer = nzb_s_inner 846 nzb_v_inner = nzb_s_inner 847 nzb_v_outer = nzb_s_inner 848 nzb_w_outer = nzb_s_inner 849 nzb_s_outer = nzb_s_inner 850 851 ! 852 !-- nzb_s_outer: 853 !-- extend nzb_local east-/westwards first, then north-/southwards 854 nzb_tmp = nzb_local 855 DO j = nys, nyn 856 DO i = nxl, nxr 857 nzb_tmp(j,i) = MAX( nzb_local(j,i-1), nzb_local(j,i), & 858 nzb_local(j,i+1) ) 714 IF ( TRIM( topography ) /= 'flat' ) THEN 715 ! 716 !-- Allocate and set the arrays containing the topography height (for output 717 !-- reasons only). 718 IF ( nxr == nx .AND. nyn /= ny ) THEN 719 ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn), & 720 zw_w_inner(nxl:nxr+1,nys:nyn) ) 721 ELSEIF ( nxr /= nx .AND. nyn == ny ) THEN 722 ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn+1), & 723 zw_w_inner(nxl:nxr,nys:nyn+1) ) 724 ELSEIF ( nxr == nx .AND. nyn == ny ) THEN 725 ALLOCATE( zu_s_inner(nxl:nxr+1,nys:nyn+1), & 726 zw_w_inner(nxl:nxr+1,nys:nyn+1) ) 727 ELSE 728 ALLOCATE( zu_s_inner(nxl:nxr,nys:nyn), & 729 zw_w_inner(nxl:nxr,nys:nyn) ) 730 ENDIF 731 732 zu_s_inner = 0.0_wp 733 zw_w_inner = 0.0_wp 734 ! 735 !-- Determine local topography height on scalar and w-grid. Note, setting 736 !-- lateral boundary values is not necessary, realized via wall_flags_static_0 737 !-- array. Further, please note that loop bounds are different from 738 !-- nxl to nxr and nys to nyn on south and right model boundary, hence, 739 !-- use intrinsic lbound and ubound functions to infer array bounds. 740 DO i = LBOUND(zu_s_inner, 1), UBOUND(zu_s_inner, 1) 741 DO j = LBOUND(zu_s_inner, 2), UBOUND(zu_s_inner, 2) 742 ! 743 !-- Topography height on scalar grid. Therefore, determine index of 744 !-- upward-facing surface element on scalar grid. 745 zu_s_inner(i,j) = zu(topo_top_ind(j,i,0)) 746 ! 747 !-- Topography height on w grid. Therefore, determine index of 748 !-- upward-facing surface element on w grid. 749 zw_w_inner(i,j) = zw(topo_top_ind(j,i,3)) 750 ENDDO 859 751 ENDDO 860 ENDDO861 862 CALL exchange_horiz_2d_int( nzb_tmp, nys, nyn, nxl, nxr, nbgp )863 864 DO i = nxl, nxr865 DO j = nys, nyn866 nzb_s_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i), &867 nzb_tmp(j+1,i) )868 ENDDO869 !870 !-- non-cyclic boundary conditions (overwritten by call of871 !-- exchange_horiz_2d_int below in case of cyclic boundary conditions)872 IF ( nys == 0 ) THEN873 j = -1874 nzb_s_outer(j,i) = MAX( nzb_tmp(j+1,i), nzb_tmp(j,i) )875 ENDIF876 IF ( nyn == ny ) THEN877 j = ny + 1878 nzb_s_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i) )879 ENDIF880 ENDDO881 !882 !-- nzb_w_outer:883 !-- identical to nzb_s_outer884 nzb_w_outer = nzb_s_outer885 !886 !-- nzb_u_inner:887 !-- extend nzb_local rightwards only888 nzb_tmp = nzb_local889 DO j = nys, nyn890 DO i = nxl, nxr891 nzb_tmp(j,i) = MAX( nzb_local(j,i-1), nzb_local(j,i) )892 ENDDO893 ENDDO894 895 CALL exchange_horiz_2d_int( nzb_tmp, nys, nyn, nxl, nxr, nbgp )896 897 nzb_u_inner = nzb_tmp898 !899 !-- nzb_u_outer:900 !-- extend current nzb_tmp (nzb_u_inner) north-/southwards901 DO i = nxl, nxr902 DO j = nys, nyn903 nzb_u_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i), &904 nzb_tmp(j+1,i) )905 ENDDO906 !907 !-- non-cyclic boundary conditions (overwritten by call of908 !-- exchange_horiz_2d_int below in case of cyclic boundary conditions)909 IF ( nys == 0 ) THEN910 j = -1911 nzb_u_outer(j,i) = MAX( nzb_tmp(j+1,i), nzb_tmp(j,i) )912 ENDIF913 IF ( nyn == ny ) THEN914 j = ny + 1915 nzb_u_outer(j,i) = MAX( nzb_tmp(j-1,i), nzb_tmp(j,i) )916 ENDIF917 ENDDO918 919 !920 !-- nzb_v_inner:921 !-- extend nzb_local northwards only922 nzb_tmp = nzb_local923 DO i = nxl, nxr924 DO j = nys, nyn925 nzb_tmp(j,i) = MAX( nzb_local(j-1,i), nzb_local(j,i) )926 ENDDO927 ENDDO928 929 CALL exchange_horiz_2d_int( nzb_tmp, nys, nyn, nxl, nxr, nbgp )930 nzb_v_inner = nzb_tmp931 932 !933 !-- nzb_v_outer:934 !-- extend current nzb_tmp (nzb_v_inner) right-/leftwards935 DO j = nys, nyn936 DO i = nxl, nxr937 nzb_v_outer(j,i) = MAX( nzb_tmp(j,i-1), nzb_tmp(j,i), &938 nzb_tmp(j,i+1) )939 ENDDO940 !941 !-- non-cyclic boundary conditions (overwritten by call of942 !-- exchange_horiz_2d_int below in case of cyclic boundary conditions)943 IF ( nxl == 0 ) THEN944 i = -1945 nzb_v_outer(j,i) = MAX( nzb_tmp(j,i+1), nzb_tmp(j,i) )946 ENDIF947 IF ( nxr == nx ) THEN948 i = nx + 1949 nzb_v_outer(j,i) = MAX( nzb_tmp(j,i-1), nzb_tmp(j,i) )950 ENDIF951 ENDDO952 953 !954 !-- Exchange of lateral boundary values (parallel computers) and cyclic955 !-- boundary conditions, if applicable.956 !-- Since nzb_s_inner and nzb_w_inner are derived directly from nzb_local957 !-- they do not require exchange and are not included here.958 CALL exchange_horiz_2d_int( nzb_u_inner, nys, nyn, nxl, nxr, nbgp )959 CALL exchange_horiz_2d_int( nzb_u_outer, nys, nyn, nxl, nxr, nbgp )960 CALL exchange_horiz_2d_int( nzb_v_inner, nys, nyn, nxl, nxr, nbgp )961 CALL exchange_horiz_2d_int( nzb_v_outer, nys, nyn, nxl, nxr, nbgp )962 CALL exchange_horiz_2d_int( nzb_w_outer, nys, nyn, nxl, nxr, nbgp )963 CALL exchange_horiz_2d_int( nzb_s_outer, nys, nyn, nxl, nxr, nbgp )964 965 !966 !-- Set the individual index arrays which define the k index from which on967 !-- the usual finite difference form (which does not use surface fluxes) is968 !-- applied969 IF ( constant_flux_layer .OR. use_surface_fluxes ) THEN970 nzb_diff_s_inner = nzb_s_inner + 2971 nzb_diff_s_outer = nzb_s_outer + 2972 ELSE973 nzb_diff_s_inner = nzb_s_inner + 1974 nzb_diff_s_outer = nzb_s_outer + 1975 752 ENDIF 976 753 ! … … 2798 2575 BTEST( wall_flags_total_0(:,:,:), ibit ) & 2799 2576 ), DIM = 1 & 2800 ) - 1 2801 2577 ) - 1 2578 2802 2579 END SUBROUTINE set_topo_flags 2803 2580 2804 2581 2805 2582 2583
Note: See TracChangeset
for help on using the changeset viewer.
