Changeset 1968 for palm/trunk/SOURCE/init_grid.f90
 Timestamp:
 Jul 18, 2016 12:01:49 PM (7 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/init_grid.f90
r1943 r1968 19 19 ! Current revisions: 20 20 !  21 ! 21 ! Changed: PEwise reading of topography file in order to avoid global definition 22 ! of arrays nzb_local and nzb_tmp. Thereby, topography definition for single 23 ! buildings and street canyons has changed, as well as flag setting for 24 ! multigrid scheme. 25 ! 26 ! Bugfix in checking l_grid anisotropy. 27 ! Simplify initial computation of lwall and vertical_influence, i.e. remove 28 ! nzb_s_inner as it is still zero at this point. 22 29 ! 23 30 ! Former revisions: … … 188 195 canyon_width_x, canyon_width_y, constant_flux_layer, & 189 196 coupling_char, dp_level_ind_b, dz, dz_max, dz_stretch_factor, & 190 dz_stretch_level, dz_stretch_level_index, ibc_uv_b, io_blocks,&191 io_ group, inflow_l, inflow_n, inflow_r, inflow_s,&197 dz_stretch_level, dz_stretch_level_index, grid_level, ibc_uv_b, & 198 io_blocks, io_group, inflow_l, inflow_n, inflow_r, inflow_s, & 192 199 masking_method, maximum_grid_level, message_string, & 193 200 momentum_advec, nest_domain, ocean, outflow_l, outflow_n, & … … 219 226 IMPLICIT NONE 220 227 221 INTEGER(iwp) :: bh !< temporary vertical index of building height 222 INTEGER(iwp) :: blx !< grid point number of building size along x 223 INTEGER(iwp) :: bly !< grid point number of building size along y 224 INTEGER(iwp) :: bxl !< index for left building wall 225 INTEGER(iwp) :: bxr !< index for right building wall 226 INTEGER(iwp) :: byn !< index for north building wall 227 INTEGER(iwp) :: bys !< index for south building wall 228 INTEGER(iwp) :: ch !< temporary vertical index for canyon height 229 INTEGER(iwp) :: cwx !< grid point number of canyon size along x 230 INTEGER(iwp) :: cwy !< grid point number of canyon size along y 231 INTEGER(iwp) :: cxl !< index for left canyon wall 232 INTEGER(iwp) :: cxr !< index for right canyon wall 233 INTEGER(iwp) :: cyn !< index for north canyon wall 234 INTEGER(iwp) :: cys !< index for south canyon wall 235 INTEGER(iwp) :: gls !< number of lateral ghost points at total model domain boundaries required for multigrid solver 236 INTEGER(iwp) :: i !< index variable along x 237 INTEGER(iwp) :: ii !< loop variable for reading topography file 238 INTEGER(iwp) :: inc !< incremental parameter for coarsening grid level 239 INTEGER(iwp) :: j !< index variable along y 240 INTEGER(iwp) :: k !< index variable along z 241 INTEGER(iwp) :: l !< loop variable 242 INTEGER(iwp) :: nxl_l !< index of left PE boundary for multigrid level 243 INTEGER(iwp) :: nxr_l !< index of right PE boundary for multigrid level 244 INTEGER(iwp) :: nyn_l !< index of north PE boundary for multigrid level 245 INTEGER(iwp) :: nys_l !< index of south PE boundary for multigrid level 246 INTEGER(iwp) :: nzb_si !< dummy index for local nzb_s_inner 247 INTEGER(iwp) :: nzt_l !< index of top PE boundary for multigrid level 248 INTEGER(iwp) :: num_hole !< number of holes (in topography) resolved by only one grid point 249 INTEGER(iwp) :: num_wall !< number of surrounding vertical walls for a single grid point 250 INTEGER(iwp) :: vi !< dummy for vertical influence 228 INTEGER(iwp) :: bh !< temporary vertical index of building height 229 INTEGER(iwp) :: blx !< grid point number of building size along x 230 INTEGER(iwp) :: bly !< grid point number of building size along y 231 INTEGER(iwp) :: bxl !< index for left building wall 232 INTEGER(iwp) :: bxr !< index for right building wall 233 INTEGER(iwp) :: byn !< index for north building wall 234 INTEGER(iwp) :: bys !< index for south building wall 235 INTEGER(iwp) :: ch !< temporary vertical index for canyon height 236 INTEGER(iwp) :: cwx !< grid point number of canyon size along x 237 INTEGER(iwp) :: cwy !< grid point number of canyon size along y 238 INTEGER(iwp) :: cxl !< index for left canyon wall 239 INTEGER(iwp) :: cxr !< index for right canyon wall 240 INTEGER(iwp) :: cyn !< index for north canyon wall 241 INTEGER(iwp) :: cys !< index for south canyon wall 242 INTEGER(iwp) :: i !< index variable along x 243 INTEGER(iwp) :: ii !< loop variable for reading topography file 244 INTEGER(iwp) :: inc !< incremental parameter for coarsening grid level 245 INTEGER(iwp) :: j !< index variable along y 246 INTEGER(iwp) :: k !< index variable along z 247 INTEGER(iwp) :: l !< loop variable 248 INTEGER(iwp) :: nxl_l !< index of left PE boundary for multigrid level 249 INTEGER(iwp) :: nxr_l !< index of right PE boundary for multigrid level 250 INTEGER(iwp) :: nyn_l !< index of north PE boundary for multigrid level 251 INTEGER(iwp) :: nys_l !< index of south PE boundary for multigrid level 252 INTEGER(iwp) :: nzb_local_max !< vertical grid index of maximum topography height 253 INTEGER(iwp) :: nzb_local_min !< vertical grid index of minimum topography height 254 INTEGER(iwp) :: nzb_si !< dummy index for local nzb_s_inner 255 INTEGER(iwp) :: nzt_l !< index of top PE boundary for multigrid level 256 INTEGER(iwp) :: num_hole !< number of holes (in topography) resolved by only one grid point 257 INTEGER(iwp) :: num_hole_l !< number of holes (in topography) resolved by only one grid point on local PE 258 INTEGER(iwp) :: num_wall !< number of surrounding vertical walls for a single grid point 259 INTEGER(iwp) :: skip_n_rows !< counting variable to skip rows while reading topography file 260 INTEGER(iwp) :: vi !< dummy for vertical influence 251 261 252 262 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: & 253 263 vertical_influence !< number of vertical grid points above obstacle where adjustment of nearwall mixing length is required 254 264 255 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: corner_nl !< index of northleft corner location to limit nearwall mixing length 256 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: corner_nr !< northright 257 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: corner_sl !< southleft 258 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: corner_sr !< southright 259 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: wall_l !< distance to adjacent leftfacing 260 !< wall 261 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: wall_n !< northfacing 262 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: wall_r !< rightfacing 263 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: wall_s !< rightfacing 264 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nzb_local !< index for topography 265 !< top at cellcenter 266 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nzb_tmp !< dummy to calculate topography indices on u and vgrid 267 268 LOGICAL :: hole = .FALSE. !< flag to check if any holes resolved by only 1 grid point were filled 269 265 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: corner_nl !< index of northleft corner location to limit nearwall mixing length 266 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: corner_nr !< northright 267 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: corner_sl !< southleft 268 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: corner_sr !< southright 269 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nzb_local !< index for topography top at cellcenter 270 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: nzb_tmp !< dummy to calculate topography indices on u and vgrid 271 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: wall_l !< distance to adjacent leftfacing wall 272 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: wall_n !< northfacing 273 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: wall_r !< rightfacing 274 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: wall_s !< rightfacing 275 276 REAL(wp) :: dum !< dummy variable to skip columns while reading topography file 270 277 REAL(wp) :: dx_l !< grid spacing along x on different multigrid level 271 278 REAL(wp) :: dy_l !< grid spacing along y on different multigrid level 272 279 REAL(wp) :: dz_stretched !< stretched vertical grid spacing 273 280 274 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: topo_height !< input variable for topography height 281 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: topo_height !< input variable for topography height 282 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zu_s_inner_l !< dummy array on global scale to write topography output array 283 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zw_w_inner_l !< dummy array on global scale to write topography output array 275 284 276 285 … … 472 481 ! Allocate outer and inner index arrays for topography and set 473 482 ! defaults. 474 ! nzb_local has to contain additional layers of ghost points for calculating475 ! the flag arrays needed for the multigrid method476 gls = 2**( maximum_grid_level )477 IF ( gls < nbgp ) gls = nbgp478 483 479 484 ALLOCATE( corner_nl(nys:nyn,nxl:nxr), corner_nr(nys:nyn,nxl:nxr), & 480 485 corner_sl(nys:nyn,nxl:nxr), corner_sr(nys:nyn,nxl:nxr), & 481 nzb_local(gls:ny+gls,gls:nx+gls), &482 nzb_tmp(nbgp:ny+nbgp,nbgp:nx+nbgp), &483 486 wall_l(nys:nyn,nxl:nxr), wall_n(nys:nyn,nxl:nxr), & 484 wall_r(nys:nyn,nxl:nxr), wall_s(nys:nyn,nxl:nxr) ) 487 wall_r(nys:nyn,nxl:nxr), wall_s(nys:nyn,nxl:nxr) ) 488 485 489 ALLOCATE( fwxm(nysg:nyng,nxlg:nxrg), fwxp(nysg:nyng,nxlg:nxrg), & 486 490 fwym(nysg:nyng,nxlg:nxrg), fwyp(nysg:nyng,nxlg:nxrg), & … … 499 503 nzb_diff_u(nysg:nyng,nxlg:nxrg), & 500 504 nzb_diff_v(nysg:nyng,nxlg:nxrg), & 505 nzb_local(nysg:nyng,nxlg:nxrg), & 506 nzb_tmp(nysg:nyng,nxlg:nxrg), & 501 507 rflags_s_inner(nzb:nzt+2,nysg:nyng,nxlg:nxrg), & 502 508 rflags_invers(nysg:nyng,nxlg:nxrg,nzb:nzt+2), & … … 559 565 ENDDO 560 566 561 DO k = 1, MAXVAL( nzb_s_inner )567 DO k = 1, nzt 562 568 IF ( l_grid(k) > 1.5_wp * dx * wall_adjustment_factor .OR. & 563 569 l_grid(k) > 1.5_wp * dy * wall_adjustment_factor ) THEN … … 573 579 vertical_influence(0) = vertical_influence(1) 574 580 575 DO i = nxlg, nxrg 576 DO j = nysg, nyng 577 DO k = nzb_s_inner(j,i) + 1, & 578 nzb_s_inner(j,i) + vertical_influence(nzb_s_inner(j,i)) 579 l_wall(k,j,i) = zu(k)  zw(nzb_s_inner(j,i)) 580 ENDDO 581 ENDDO 581 DO k = nzb + 1, nzb + vertical_influence(nzb) 582 l_wall(k,:,:) = zu(k)  zw(nzb) 582 583 ENDDO 583 584 … … 629 630 630 631 ! 631 ! Define the building. 632 ! Define the building. 632 633 nzb_local = 0 633 nzb_local(bys:byn,bxl:bxr) = bh 634 IF ( bxl <= nxr .AND. bxr >= nxl .AND. & 635 bys <= nyn .AND. byn >= nys ) & 636 nzb_local(MAX(nys,bys):MIN(nyn,byn),MAX(nxl,bxl):MIN(nxr,bxr)) = bh 634 637 635 638 CASE ( 'single_street_canyon' ) … … 700 703 nzb_local = ch 701 704 IF ( canyon_width_x /= 9999999.9_wp ) THEN 702 nzb_local(:,cxl+1:cxr1) = 0 705 IF ( cxl <= nxr .AND. cxr >= nxl ) & 706 nzb_local(:,MAX(nxl,cxl+1):MIN(nxr,cxr1)) = 0 703 707 ELSEIF ( canyon_width_y /= 9999999.9_wp ) THEN 704 nzb_local(cys+1:cyn1,:) = 0 708 IF ( cys <= nyn .AND. cyn >= nys ) & 709 nzb_local(MAX(nys,cys+1):MIN(nyn,cyn1),:) = 0 705 710 ENDIF 706 711 707 712 CASE ( 'read_from_file' ) 708 713 709 ALLOCATE ( topo_height( 0:ny,0:nx) )714 ALLOCATE ( topo_height(nys:nyn,nxl:nxr) ) 710 715 711 716 DO ii = 0, io_blocks1 … … 717 722 OPEN( 90, FILE='TOPOGRAPHY_DATA'//TRIM( coupling_char ), & 718 723 STATUS='OLD', FORM='FORMATTED', ERR=10 ) 719 DO j = ny, 0, 1 720 READ( 90, *, ERR=11, END=11 ) ( topo_height(j,i), i = 0,nx ) 724 ! 725 ! Read topography PEwise. Rows are read from nyn to nys, columns 726 ! are read from nxl to nxr. At first, nynyn rows need to be skipped. 727 skip_n_rows = 0 728 DO WHILE ( skip_n_rows < ny  nyn ) 729 READ( 90, * ) 730 skip_n_rows = skip_n_rows + 1 731 ENDDO 732 ! 733 ! Read data from nyn to nys and nxl to nxr. Therefore, skip 734 ! column until nxl1 is reached 735 DO j = nyn, nys, 1 736 READ( 90, *, ERR=11, END=11 ) & 737 ( dum, i = 0, nxl1 ), & 738 ( topo_height(j,i), i = nxl, nxr ) 721 739 ENDDO 722 740 … … 741 759 ! 742 760 ! Calculate the index height of the topography 743 DO i = 0, nx 744 DO j = 0, ny 761 nzb_local = 0 762 DO i = nxl, nxr 763 DO j = nys, nyn 745 764 nzb_local(j,i) = MINLOC( ABS( zw  topo_height(j,i) ), 1 )  1 746 765 IF ( ABS( zw(nzb_local(j,i) )  topo_height(j,i) ) == & … … 759 778 ! Before checking for holes, set lateral boundary conditions for 760 779 ! topography. After holefilling, boundary conditions must be set again! 761 IF ( bc_ns_cyc ) THEN 762 nzb_local(1,:) = nzb_local(ny,:) 763 nzb_local(ny+1,:) = nzb_local(0,:) 764 ELSE 765 nzb_local(1,:) = nzb_local(0,:) 766 nzb_local(ny+1,:) = nzb_local(ny,:) 767 ENDIF 768 769 IF ( bc_lr_cyc ) THEN 770 nzb_local(:,1) = nzb_local(:,nx) 771 nzb_local(:,nx+1) = nzb_local(:,0) 772 ELSE 773 nzb_local(:,1) = nzb_local(:,0) 774 nzb_local(:,nx+1) = nzb_local(:,nx) 775 ENDIF 776 777 num_hole = 0 778 DO i = 0, nx 779 DO j = 0, ny 780 CALL exchange_horiz_2d_int( nzb_local, nys, nyn, nxl, nxr, nbgp ) 781 782 IF ( .NOT. bc_ns_cyc ) THEN 783 IF ( nys == 0 ) nzb_local(1,:) = nzb_local(0,:) 784 IF ( nyn == ny ) nzb_local(ny+1,:) = nzb_local(ny,:) 785 ENDIF 786 787 IF ( .NOT. bc_lr_cyc ) THEN 788 IF ( nxl == 0 ) nzb_local(:,1) = nzb_local(:,0) 789 IF ( nxr == nx ) nzb_local(:,nx+1) = nzb_local(:,nx) 790 ENDIF 791 792 num_hole_l = 0 793 DO i = nxl, nxr 794 DO j = nys, nyn 780 795 781 796 num_wall = 0 … … 791 806 792 807 IF ( num_wall == 4 ) THEN 793 hole = .TRUE.794 808 nzb_local(j,i) = MIN( nzb_local(j1,i), nzb_local(j+1,i), & 795 809 nzb_local(j,i1), nzb_local(j,i+1) ) 796 num_hole = num_hole+ 1810 num_hole_l = num_hole_l + 1 797 811 ENDIF 798 812 ENDDO 799 813 ENDDO 800 814 ! 815 ! Count the total number of holes, required for informative message. 816 #if defined( __parallel ) 817 CALL MPI_ALLREDUCE( num_hole_l, num_hole, 1, MPI_INTEGER, MPI_SUM, & 818 comm2d, ierr ) 819 #else 820 num_hole = num_hole_l 821 #endif 822 ! 801 823 ! Create an informative message if any hole was removed. 802 IF ( hole) THEN824 IF ( num_hole > 0 ) THEN 803 825 WRITE( message_string, * ) num_hole, 'hole(s) resolved by only '//& 804 826 'one grid point were filled' … … 806 828 ENDIF 807 829 ! 808 ! Add cyclic or Neumann boundary conditions (additional layers are for 809 ! calculating flag arrays needed for the multigrid sover) 810 IF ( bc_ns_cyc ) THEN 811 nzb_local(gls:1,0:nx) = nzb_local(nygls+1:ny,0:nx) 812 nzb_local(ny+1:ny+gls,0:nx) = nzb_local(0:gls1,0:nx) 813 ELSE 814 DO j = gls, 1 815 nzb_local(j,0:nx) = nzb_local(0,0:nx) 816 ENDDO 817 DO j = ny+1, ny+gls 818 nzb_local(j,0:nx) = nzb_local(ny,0:nx) 819 ENDDO 820 ENDIF 821 822 IF ( bc_lr_cyc ) THEN 823 nzb_local(:,gls:1) = nzb_local(:,nxgls+1:nx) 824 nzb_local(:,nx+1:nx+gls) = nzb_local(:,0:gls1) 825 ELSE 826 DO i = gls, 1 827 nzb_local(:,i) = nzb_local(:,0) 828 ENDDO 829 DO i = nx+1, nx+gls 830 nzb_local(:,i) = nzb_local(:,nx) 831 ENDDO 830 ! Exchange ghostpoints, as well as add cyclic or Neumann boundary 831 ! conditions. 832 CALL exchange_horiz_2d_int( nzb_local, nys, nyn, nxl, nxr, nbgp ) 833 834 IF ( .NOT. bc_ns_cyc ) THEN 835 IF ( nys == 0 ) nzb_local(1,:) = nzb_local(0,:) 836 IF ( nyn == ny ) nzb_local(ny+1,:) = nzb_local(ny,:) 837 ENDIF 838 839 IF ( .NOT. bc_lr_cyc ) THEN 840 IF ( nxl == 0 ) nzb_local(:,1) = nzb_local(:,0) 841 IF ( nxr == nx ) nzb_local(:,nx+1) = nzb_local(:,nx) 832 842 ENDIF 833 843 … … 838 848 ! case in the user interface. There, the subroutine user_init_grid 839 849 ! checks which of these two conditions applies. 840 CALL user_init_grid( gls,nzb_local )850 CALL user_init_grid( nzb_local ) 841 851 842 852 END SELECT … … 845 855 ! steering the degradation of order of the applied advection scheme. 846 856 ! In case of noncyclic lateral boundaries, the order of the advection 847 ! scheme have to be reduced up to nzt (required at the lateral boundaries). 857 ! scheme has to be reduced up to nzt (required at the lateral boundaries). 858 #if defined( __parallel ) 859 CALL MPI_ALLREDUCE( MAXVAL( nzb_local ) + 1, nzb_max, 1, MPI_INTEGER, & 860 MPI_MAX, comm2d, ierr ) 861 #else 848 862 nzb_max = MAXVAL( nzb_local ) + 1 863 #endif 849 864 IF ( inflow_l .OR. outflow_l .OR. inflow_r .OR. outflow_r .OR. & 850 865 inflow_n .OR. outflow_n .OR. inflow_s .OR. outflow_s .OR. & … … 859 874 ! zu_s_inner and zw_w_inner 860 875 IF ( TRIM( topography ) /= 'flat' ) THEN 861 876 #if defined( __parallel ) 877 CALL MPI_ALLREDUCE( MAXVAL( nzb_local ), nzb_local_max, 1, MPI_INTEGER, & 878 MPI_MAX, comm2d, ierr ) 879 CALL MPI_ALLREDUCE( MAXVAL( nzb_local ), nzb_local_min, 1, MPI_INTEGER, & 880 MPI_MIN, comm2d, ierr ) 881 #else 882 nzb_local_max = MAXVAL( nzb_local ) 883 nzb_local_min = MINVAL( nzb_local ) 884 #endif 862 885 ! 863 886 ! Consistency checks 864 IF ( MINVAL( nzb_local ) < 0 .OR. MAXVAL( nzb_local )> nz + 1 ) THEN887 IF ( nzb_local_min < 0 .OR. nzb_local_max > nz + 1 ) THEN 865 888 WRITE( message_string, * ) 'nzb_local values are outside the', & 866 889 'model domain', & 867 '&MINVAL( nzb_local ) = ', MINVAL(nzb_local),&868 '&MAXVAL( nzb_local ) = ', MAXVAL(nzb_local)890 '&MINVAL( nzb_local ) = ', nzb_local_min, & 891 '&MAXVAL( nzb_local ) = ', nzb_local_max 869 892 CALL message( 'init_grid', 'PA0210', 1, 2, 0, 6, 0 ) 870 ENDIF871 872 IF ( bc_lr_cyc ) THEN873 IF ( ANY( nzb_local(:,1) /= nzb_local(:,nx) ) .OR. &874 ANY( nzb_local(:,0) /= nzb_local(:,nx+1) ) ) THEN875 message_string = 'nzb_local does not fulfill cyclic' // &876 ' boundary condition in xdirection'877 CALL message( 'init_grid', 'PA0211', 1, 2, 0, 6, 0 )878 ENDIF879 ENDIF880 IF ( bc_ns_cyc ) THEN881 IF ( ANY( nzb_local(1,:) /= nzb_local(ny,:) ) .OR. &882 ANY( nzb_local(0,:) /= nzb_local(ny+1,:) ) ) THEN883 message_string = 'nzb_local does not fulfill cyclic' // &884 ' boundary condition in ydirection'885 CALL message( 'init_grid', 'PA0212', 1, 2, 0, 6, 0 )886 ENDIF887 893 ENDIF 888 894 … … 896 902 ! Therefore, the extent of topography in nzb_local is now reduced by 897 903 ! 1dx at the E topography walls and by 1dy at the N topography walls 898 ! to form the basis for nzb_s_inner. 899 DO j = gls, ny + gls 900 DO i = gls, nx 904 ! to form the basis for nzb_s_inner. 905 ! Note, the reverse memory access (ij instead of ji) is absolutely 906 ! required at this point. 907 DO j = nys+1, nyn+1 908 DO i = nxl1, nxr 901 909 nzb_local(j,i) = MIN( nzb_local(j,i), nzb_local(j,i+1) ) 902 910 ENDDO 903 911 ENDDO 904 ! apply cyclic boundary conditions in xdirection 905 !(ist das erforderlich? Ursache von Seung Bus Fehler?) 906 nzb_local(:,nx+1:nx+gls) = nzb_local(:,0:gls1) 907 DO i = gls, nx + gls 908 DO j = gls, ny 912 ! 913 ! Exchange ghost points 914 CALL exchange_horiz_2d_int( nzb_local, nys, nyn, nxl, nxr, nbgp ) 915 916 DO i = nxl, nxr+1 917 DO j = nys1, nyn 909 918 nzb_local(j,i) = MIN( nzb_local(j,i), nzb_local(j+1,i) ) 910 919 ENDDO 911 920 ENDDO 912 !  apply cyclic boundary conditions in ydirection913 ! (ist das erforderlich? Ursache von Seung Bus Fehler?)914 nzb_local(ny+1:ny+gls,:) = nzb_local(0:gls1,:)921 ! 922 ! Exchange ghost points 923 CALL exchange_horiz_2d_int( nzb_local, nys, nyn, nxl, nxr, nbgp ) 915 924 ENDIF 916 917 925 ! 918 926 ! Initialize index arrays nzb_s_inner and nzb_w_inner 919 nzb_s_inner = nzb_local (nysg:nyng,nxlg:nxrg)920 nzb_w_inner = nzb_local (nysg:nyng,nxlg:nxrg)927 nzb_s_inner = nzb_local 928 nzb_w_inner = nzb_local 921 929 922 930 ! … … 937 945 ! nzb_s_outer: 938 946 ! extend nzb_local east/westwards first, then north/southwards 939 nzb_tmp = nzb_local (nbgp:ny+nbgp,nbgp:nx+nbgp)940 DO j = 1, ny + 1941 DO i = 0, nx947 nzb_tmp = nzb_local 948 DO j = nys, nyn 949 DO i = nxl, nxr 942 950 nzb_tmp(j,i) = MAX( nzb_local(j,i1), nzb_local(j,i), & 943 951 nzb_local(j,i+1) ) 944 952 ENDDO 945 953 ENDDO 954 955 CALL exchange_horiz_2d_int( nzb_tmp, nys, nyn, nxl, nxr, nbgp ) 956 946 957 DO i = nxl, nxr 947 958 DO j = nys, nyn … … 969 980 ! nzb_u_inner: 970 981 ! extend nzb_local rightwards only 971 nzb_tmp = nzb_local (nbgp:ny+nbgp,nbgp:nx+nbgp)972 DO j = 1, ny + 1973 DO i = 0, nx + 1982 nzb_tmp = nzb_local 983 DO j = nys, nyn 984 DO i = nxl, nxr 974 985 nzb_tmp(j,i) = MAX( nzb_local(j,i1), nzb_local(j,i) ) 975 986 ENDDO 976 987 ENDDO 977 nzb_u_inner = nzb_tmp(nysg:nyng,nxlg:nxrg) 978 988 989 CALL exchange_horiz_2d_int( nzb_tmp, nys, nyn, nxl, nxr, nbgp ) 990 991 nzb_u_inner = nzb_tmp 979 992 ! 980 993 ! nzb_u_outer: … … 1001 1014 ! nzb_v_inner: 1002 1015 ! extend nzb_local northwards only 1003 nzb_tmp = nzb_local (nbgp:ny+nbgp,nbgp:nx+nbgp)1004 DO i = 1, nx + 11005 DO j = 0, ny + 11016 nzb_tmp = nzb_local 1017 DO i = nxl, nxr 1018 DO j = nys, nyn 1006 1019 nzb_tmp(j,i) = MAX( nzb_local(j1,i), nzb_local(j,i) ) 1007 1020 ENDDO 1008 1021 ENDDO 1009 nzb_v_inner = nzb_tmp(nysg:nyng,nxlg:nxrg) 1022 1023 CALL exchange_horiz_2d_int( nzb_tmp, nys, nyn, nxl, nxr, nbgp ) 1024 nzb_v_inner = nzb_tmp 1010 1025 1011 1026 ! … … 1035 1050 ! Since nzb_s_inner and nzb_w_inner are derived directly from nzb_local 1036 1051 ! they do not require exchange and are not included here. 1037 CALL exchange_horiz_2d_int( nzb_u_inner )1038 CALL exchange_horiz_2d_int( nzb_u_outer )1039 CALL exchange_horiz_2d_int( nzb_v_inner )1040 CALL exchange_horiz_2d_int( nzb_v_outer )1041 CALL exchange_horiz_2d_int( nzb_w_outer )1042 CALL exchange_horiz_2d_int( nzb_s_outer )1052 CALL exchange_horiz_2d_int( nzb_u_inner, nys, nyn, nxl, nxr, nbgp ) 1053 CALL exchange_horiz_2d_int( nzb_u_outer, nys, nyn, nxl, nxr, nbgp ) 1054 CALL exchange_horiz_2d_int( nzb_v_inner, nys, nyn, nxl, nxr, nbgp ) 1055 CALL exchange_horiz_2d_int( nzb_v_outer, nys, nyn, nxl, nxr, nbgp ) 1056 CALL exchange_horiz_2d_int( nzb_w_outer, nys, nyn, nxl, nxr, nbgp ) 1057 CALL exchange_horiz_2d_int( nzb_s_outer, nys, nyn, nxl, nxr, nbgp ) 1043 1058 1044 1059 ! 1045 1060 ! Allocate and set the arrays containing the topography height 1046 IF ( myid == 0 ) THEN 1047 1048 ALLOCATE( zu_s_inner(0:nx+1,0:ny+1), zw_w_inner(0:nx+1,0:ny+1) ) 1049 1050 DO i = 0, nx + 1 1051 DO j = 0, ny + 1 1052 zu_s_inner(i,j) = zu(nzb_local(j,i)) 1053 zw_w_inner(i,j) = zw(nzb_local(j,i)) 1054 ENDDO 1061 ALLOCATE( zu_s_inner(0:nx+1,0:ny+1), zw_w_inner(0:nx+1,0:ny+1), & 1062 zu_s_inner_l(0:nx+1,0:ny+1), zw_w_inner_l(0:nx+1,0:ny+1) ) 1063 1064 zu_s_inner = 0.0_wp 1065 zw_w_inner = 0.0_wp 1066 zu_s_inner_l = 0.0_wp 1067 zw_w_inner_l = 0.0_wp 1068 1069 DO i = nxl, nxr 1070 DO j = nys, nyn 1071 zu_s_inner_l(i,j) = zu(nzb_local(j,i)) 1072 zw_w_inner_l(i,j) = zw(nzb_local(j,i)) 1055 1073 ENDDO 1056 1057 ENDIF 1074 ENDDO 1075 1076 #if defined( __parallel ) 1077 CALL MPI_REDUCE( zu_s_inner_l, zu_s_inner, (nx+2)*(ny+2), & 1078 MPI_REAL, MPI_SUM, 0, comm2d, ierr ) 1079 CALL MPI_REDUCE( zw_w_inner_l, zw_w_inner, (nx+2)*(ny+2), & 1080 MPI_REAL, MPI_SUM, 0, comm2d, ierr ) 1081 #else 1082 zu_s_inner = zu_s_inner_l 1083 zw_w_inner = zw_w_inner_l 1084 #endif 1085 1086 DEALLOCATE( zu_s_inner_l, zw_w_inner_l ) 1087 IF ( myid /= 0 ) DEALLOCATE( zu_s_inner, zw_w_inner ) 1088 ! 1089 ! Set south and left ghost points, required for netcdf output 1090 IF ( myid == 0 ) THEN 1091 IF( bc_lr_cyc ) THEN 1092 zu_s_inner(nx+1,:) = zu_s_inner(0,:) 1093 zw_w_inner(nx+1,:) = zw_w_inner(0,:) 1094 ELSE 1095 zu_s_inner(nx+1,:) = zu_s_inner(nx,:) 1096 zw_w_inner(nx+1,:) = zw_w_inner(nx,:) 1097 ENDIF 1098 IF( bc_ns_cyc ) THEN 1099 zu_s_inner(:,ny+1) = zu_s_inner(:,0) 1100 zw_w_inner(:,ny+1) = zw_w_inner(:,0) 1101 ELSE 1102 zu_s_inner(:,ny+1) = zu_s_inner(:,ny) 1103 zw_w_inner(:,ny+1) = zw_w_inner(:,ny) 1104 ENDIF 1105 ENDIF 1058 1106 ! 1059 1107 ! Set flag arrays to be used for masking of grid points … … 1068 1116 1069 1117 ENDIF 1118 ! 1119 ! Deallocate temporary array, as it might be reused for different 1120 ! gridlevels further below. 1121 DEALLOCATE( nzb_tmp ) 1070 1122 1071 1123 ! … … 1190 1242 ENDDO 1191 1243 ENDDO 1192 1193 1244 ! 1194 1245 ! Calculate wall flag arrays for the multigrid method. … … 1196 1247 ! version. 1197 1248 IF ( psolver == 'multigrid_noopt' ) THEN 1198 ! 1199 ! Gridpoint increment of the current level 1249 1250 ! 1251 ! Gridpoint increment of the current level. 1200 1252 inc = 1 1201 1202 1253 DO l = maximum_grid_level, 1 , 1 1254 ! 1255 ! Set grid_level as it is required for exchange_horiz_2d_int 1256 grid_level = l 1203 1257 1204 1258 nxl_l = nxl_mg(l) … … 1207 1261 nyn_l = nyn_mg(l) 1208 1262 nzt_l = nzt_mg(l) 1209 1210 1263 ! 1211 1264 ! Assign the flag level to be calculated … … 1251 1304 IF ( .NOT. masking_method ) THEN 1252 1305 1306 ! 1307 ! Allocate temporary array for topography heights on coarser grid 1308 ! level. Please note, 2 ghoist points are required, in order to 1309 ! calculate flags() on the interior ghost point. 1310 ALLOCATE( nzb_tmp(nys_l2:nyn_l+2,nxl_l2:nxr_l+2) ) 1311 nzb_tmp = 0 1312 1313 DO i = nxl_l, nxr_l 1314 DO j = nys_l, nyn_l 1315 nzb_tmp(j,i) = nzb_local(j*inc,i*inc) 1316 ENDDO 1317 ENDDO 1318 ! 1319 ! Exchange ghost points on respective multigrid level. 2 ghost points 1320 ! are required, in order to calculate flags on 1321 ! nys_l1 / nyn_l+1 / nxl_l1 / nxr_l+1. The alternative would be to 1322 ! exchange 3DINTEGER array flags on the respective multigrid level. 1323 CALL exchange_horiz_2d_int( nzb_tmp, nys_l, nyn_l, nxl_l, nxr_l, 2 ) 1324 ! 1325 ! Set noncyclic boundary conditions on respective multigrid level 1326 IF ( .NOT. bc_ns_cyc ) THEN 1327 IF ( nys == 0 ) THEN 1328 nzb_tmp(2,:) = nzb_tmp(0,:) 1329 nzb_tmp(1,:) = nzb_tmp(0,:) 1330 ENDIF 1331 IF ( nyn == ny ) THEN 1332 nzb_tmp(ny+2,:) = nzb_tmp(ny,:) 1333 nzb_tmp(ny+1,:) = nzb_tmp(ny,:) 1334 ENDIF 1335 ENDIF 1336 IF ( .NOT. bc_lr_cyc ) THEN 1337 IF ( nxl == 0 ) THEN 1338 nzb_tmp(:,2) = nzb_tmp(:,0) 1339 nzb_tmp(:,1) = nzb_tmp(:,0) 1340 ENDIF 1341 IF ( nxr == nx ) THEN 1342 nzb_tmp(:,nx+1) = nzb_tmp(:,nx) 1343 nzb_tmp(:,nx+2) = nzb_tmp(:,nx) 1344 ENDIF 1345 ENDIF 1346 1253 1347 DO i = nxl_l1, nxr_l+1 1254 1348 DO j = nys_l1, nyn_l+1 1255 DO k = nzb, nzt_l+1 1256 1349 DO k = nzb, nzt_l+1 1257 1350 ! 1258 1351 ! Inside/outside building (inside building does not need 1259 1352 ! further tests for walls) 1260 IF ( k*inc <= nzb_ local(j*inc,i*inc) ) THEN1353 IF ( k*inc <= nzb_tmp(j,i) ) THEN 1261 1354 1262 1355 flags(k,j,i) = IBSET( flags(k,j,i), 6 ) … … 1265 1358 ! 1266 1359 ! Bottom wall 1267 IF ( (k1)*inc <= nzb_ local(j*inc,i*inc) ) THEN1360 IF ( (k1)*inc <= nzb_tmp(j,i) ) THEN 1268 1361 flags(k,j,i) = IBSET( flags(k,j,i), 0 ) 1269 1362 ENDIF 1270 1363 ! 1271 1364 ! South wall 1272 IF ( k*inc <= nzb_ local((j1)*inc,i*inc) ) THEN1365 IF ( k*inc <= nzb_tmp(j1,i) ) THEN 1273 1366 flags(k,j,i) = IBSET( flags(k,j,i), 2 ) 1274 1367 ENDIF 1275 1368 ! 1276 1369 ! North wall 1277 IF ( k*inc <= nzb_ local((j+1)*inc,i*inc) ) THEN1370 IF ( k*inc <= nzb_tmp(j+1,i) ) THEN 1278 1371 flags(k,j,i) = IBSET( flags(k,j,i), 3 ) 1279 1372 ENDIF 1280 1373 ! 1281 1374 ! Left wall 1282 IF ( k*inc <= nzb_ local(j*inc,(i1)*inc) ) THEN1375 IF ( k*inc <= nzb_tmp(j,i1) ) THEN 1283 1376 flags(k,j,i) = IBSET( flags(k,j,i), 4 ) 1284 1377 ENDIF 1285 1378 ! 1286 1379 ! Right wall 1287 IF ( k*inc <= nzb_ local(j*inc,(i+1)*inc) ) THEN1380 IF ( k*inc <= nzb_tmp(j,i+1) ) THEN 1288 1381 flags(k,j,i) = IBSET( flags(k,j,i), 5 ) 1289 1382 ENDIF … … 1295 1388 ENDDO 1296 1389 1390 DEALLOCATE( nzb_tmp ) 1391 1297 1392 ENDIF 1298 1393 … … 1300 1395 1301 1396 ENDDO 1302 1397 ! 1398 ! Reset grid_level to "normal" grid 1399 grid_level = 0 1400 1303 1401 ENDIF 1304 1402 ! … … 1464 1562 1465 1563 DEALLOCATE( corner_nl, corner_nr, corner_sl, corner_sr, nzb_local, & 1466 nzb_tmp,vertical_influence, wall_l, wall_n, wall_r, wall_s )1564 vertical_influence, wall_l, wall_n, wall_r, wall_s ) 1467 1565 1468 1566
Note: See TracChangeset
for help on using the changeset viewer.