Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Ignore:
Timestamp:
Jul 18, 2016 12:01:49 PM (9 years ago)
Author:
suehring
Message:

PE-wise reading of topography file

File:
1 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/init_grid.f90

    r1943 r1968  
    1919! Current revisions:
    2020! -----------------
    21 !
     21! Changed: PE-wise 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.
    2229!
    2330! Former revisions:
    … …  
    188195               canyon_width_x, canyon_width_y, constant_flux_layer,            &
    189196               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,    &
    192199               masking_method, maximum_grid_level, message_string,             &
    193200               momentum_advec, nest_domain, ocean, outflow_l, outflow_n,       &
    … …  
    219226    IMPLICIT NONE
    220227
    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
    251261
    252262    INTEGER(iwp), DIMENSION(:), ALLOCATABLE   ::                               &
    253263                     vertical_influence  !< number of vertical grid points above obstacle where adjustment of near-wall mixing length is required
    254264                                         
    255     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  corner_nl  !< index of north-left corner location to limit near-wall mixing length
    256     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  corner_nr  !< north-right
    257     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  corner_sl  !< south-left
    258     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  corner_sr  !< south-right
    259     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  wall_l     !< distance to adjacent left-facing
    260                                                              !< wall
    261     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  wall_n     !< north-facing
    262     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  wall_r     !< right-facing
    263     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  wall_s     !< right-facing
    264     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_local  !< index for topography
    265                                                              !< top at cell-center
    266     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_tmp    !< dummy to calculate topography indices on u- and v-grid
    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 north-left corner location to limit near-wall mixing length
     266    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  corner_nr      !< north-right
     267    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  corner_sl      !< south-left
     268    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  corner_sr      !< south-right
     269    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_local      !< index for topography top at cell-center
     270    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  nzb_tmp        !< dummy to calculate topography indices on u- and v-grid
     271    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  wall_l         !< distance to adjacent left-facing wall
     272    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  wall_n         !< north-facing
     273    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  wall_r         !< right-facing
     274    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  wall_s         !< right-facing
     275
     276    REAL(wp) ::  dum           !< dummy variable to skip columns while reading topography file   
    270277    REAL(wp) ::  dx_l          !< grid spacing along x on different multigrid level
    271278    REAL(wp) ::  dy_l          !< grid spacing along y on different multigrid level
    272279    REAL(wp) ::  dz_stretched  !< stretched vertical grid spacing
    273280
    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
    275284
    276285   
    … …  
    472481!-- Allocate outer and inner index arrays for topography and set
    473482!-- defaults.
    474 !-- nzb_local has to contain additional layers of ghost points for calculating
    475 !-- the flag arrays needed for the multigrid method
    476     gls = 2**( maximum_grid_level )
    477     IF ( gls < nbgp )  gls = nbgp
    478483
    479484    ALLOCATE( corner_nl(nys:nyn,nxl:nxr), corner_nr(nys:nyn,nxl:nxr),       &
    480485              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),                         &
    483486              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     
    485489    ALLOCATE( fwxm(nysg:nyng,nxlg:nxrg), fwxp(nysg:nyng,nxlg:nxrg),         &
    486490              fwym(nysg:nyng,nxlg:nxrg), fwyp(nysg:nyng,nxlg:nxrg),         &
    … …  
    499503              nzb_diff_u(nysg:nyng,nxlg:nxrg),                              &
    500504              nzb_diff_v(nysg:nyng,nxlg:nxrg),                              &
     505              nzb_local(nysg:nyng,nxlg:nxrg),                               &
     506              nzb_tmp(nysg:nyng,nxlg:nxrg),                                 &
    501507              rflags_s_inner(nzb:nzt+2,nysg:nyng,nxlg:nxrg),                &
    502508              rflags_invers(nysg:nyng,nxlg:nxrg,nzb:nzt+2),                 &
    … …  
    559565    ENDDO
    560566
    561     DO  k = 1, MAXVAL( nzb_s_inner )
     567    DO  k = 1, nzt
    562568       IF ( l_grid(k) > 1.5_wp * dx * wall_adjustment_factor .OR.  &
    563569            l_grid(k) > 1.5_wp * dy * wall_adjustment_factor )  THEN
    … …  
    573579    vertical_influence(0) = vertical_influence(1)
    574580
    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)
    582583    ENDDO
    583584
    … …  
    629630
    630631!
    631 !--       Define the building.
     632!--       Define the building. 
    632633          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
    634637
    635638       CASE ( 'single_street_canyon' )
    … …  
    700703          nzb_local = ch
    701704          IF ( canyon_width_x /= 9999999.9_wp )  THEN
    702              nzb_local(:,cxl+1:cxr-1) = 0
     705             IF ( cxl <= nxr  .AND.  cxr >= nxl )                              &
     706                nzb_local(:,MAX(nxl,cxl+1):MIN(nxr,cxr-1)) = 0
    703707          ELSEIF ( canyon_width_y /= 9999999.9_wp )  THEN
    704              nzb_local(cys+1:cyn-1,:) = 0
     708             IF ( cys <= nyn  .AND.  cyn >= nys )                              &         
     709                nzb_local(MAX(nys,cys+1):MIN(nyn,cyn-1),:) = 0
    705710          ENDIF
    706711
    707712       CASE ( 'read_from_file' )
    708713
    709           ALLOCATE ( topo_height(0:ny,0:nx) )
     714          ALLOCATE ( topo_height(nys:nyn,nxl:nxr) )
    710715
    711716          DO  ii = 0, io_blocks-1
    … …  
    717722                OPEN( 90, FILE='TOPOGRAPHY_DATA'//TRIM( coupling_char ),       &
    718723                          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 PE-wise. Rows are read from nyn to nys, columns
     726!--             are read from nxl to nxr. At first, ny-nyn 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 nxl-1 is reached
     735                DO  j = nyn, nys, -1
     736                   READ( 90, *, ERR=11, END=11 )                               &
     737                                              ( dum, i = 0, nxl-1 ),           &
     738                                              ( topo_height(j,i), i = nxl, nxr )
    721739                ENDDO
    722740
    … …  
    741759!
    742760!--       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
    745764                nzb_local(j,i) = MINLOC( ABS( zw - topo_height(j,i) ), 1 ) - 1
    746765                IF ( ABS( zw(nzb_local(j,i)  ) - topo_height(j,i) ) == &
    … …  
    759778!--       Before checking for holes, set lateral boundary conditions for
    760779!--       topography. After hole-filling, 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
    780795
    781796                num_wall = 0
    … …  
    791806
    792807                IF ( num_wall == 4 )  THEN
    793                    hole           = .TRUE.
    794808                   nzb_local(j,i) = MIN( nzb_local(j-1,i), nzb_local(j+1,i),   &
    795809                                         nzb_local(j,i-1), nzb_local(j,i+1) )
    796                    num_hole       = num_hole + 1
     810                   num_hole_l     = num_hole_l + 1
    797811                ENDIF
    798812             ENDDO
    799813          ENDDO
    800814!
     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!
    801823!--       Create an informative message if any hole was removed.
    802           IF ( hole )  THEN
     824          IF ( num_hole > 0 )  THEN
    803825             WRITE( message_string, * ) num_hole, 'hole(s) resolved by only '//&
    804826                                                  'one grid point were filled'
    … …  
    806828          ENDIF
    807829!
    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(ny-gls+1:ny,0:nx)
    812              nzb_local(ny+1:ny+gls,0:nx) = nzb_local(0:gls-1,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(:,nx-gls+1:nx)
    824              nzb_local(:,nx+1:nx+gls) = nzb_local(:,0:gls-1)
    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 ghost-points, 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)         
    832842          ENDIF
    833843
    … …  
    838848!--       case in the user interface. There, the subroutine user_init_grid
    839849!--       checks which of these two conditions applies.
    840           CALL user_init_grid( gls, nzb_local )
     850          CALL user_init_grid( nzb_local )
    841851
    842852    END SELECT
    … …  
    845855!-- steering the degradation of order of the applied advection scheme.
    846856!-- In case of non-cyclic 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
    848862    nzb_max = MAXVAL( nzb_local ) + 1
     863#endif
    849864    IF ( inflow_l .OR. outflow_l .OR. inflow_r .OR. outflow_r .OR.             &
    850865         inflow_n .OR. outflow_n .OR. inflow_s .OR. outflow_s .OR.             &
    … …  
    859874!-- zu_s_inner and zw_w_inner
    860875    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
    862885!
    863886!--    Consistency checks
    864        IF ( MINVAL( nzb_local ) < 0  .OR.  MAXVAL( nzb_local ) > nz + 1 )  THEN
     887       IF ( nzb_local_min < 0  .OR.  nzb_local_max > nz + 1 )  THEN
    865888          WRITE( message_string, * ) 'nzb_local values are outside the',       &
    866889                                '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
    869892          CALL message( 'init_grid', 'PA0210', 1, 2, 0, 6, 0 )
    870        ENDIF
    871 
    872        IF ( bc_lr_cyc )  THEN
    873           IF ( ANY( nzb_local(:,-1) /= nzb_local(:,nx)   )  .OR.               &
    874                ANY( nzb_local(:,0)  /= nzb_local(:,nx+1) ) )  THEN
    875              message_string = 'nzb_local does not fulfill cyclic' //           &
    876                               ' boundary condition in x-direction'
    877              CALL message( 'init_grid', 'PA0211', 1, 2, 0, 6, 0 )
    878           ENDIF
    879        ENDIF
    880        IF ( bc_ns_cyc )  THEN
    881           IF ( ANY( nzb_local(-1,:) /= nzb_local(ny,:)   )  .OR.               &
    882                ANY( nzb_local(0,:)  /= nzb_local(ny+1,:) ) )  THEN
    883              message_string = 'nzb_local does not fulfill cyclic' //           &
    884                               ' boundary condition in y-direction'
    885              CALL message( 'init_grid', 'PA0212', 1, 2, 0, 6, 0 )
    886           ENDIF
    887893       ENDIF
    888894
    … …  
    896902!--       Therefore, the extent of topography in nzb_local is now reduced by
    897903!--       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 (i-j instead of j-i) is absolutely
     906!--       required at this point.
     907          DO  j = nys+1, nyn+1
     908             DO  i = nxl-1, nxr
    901909                nzb_local(j,i) = MIN( nzb_local(j,i), nzb_local(j,i+1) )
    902910             ENDDO
    903911          ENDDO
    904 !--       apply cyclic boundary conditions in x-direction
    905 !(ist das erforderlich? Ursache von Seung Bus Fehler?)
    906           nzb_local(:,nx+1:nx+gls) = nzb_local(:,0:gls-1)
    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 = nys-1, nyn
    909918                nzb_local(j,i) = MIN( nzb_local(j,i), nzb_local(j+1,i) )
    910919             ENDDO
    911920          ENDDO
    912 !--       apply cyclic boundary conditions in y-direction
    913 !(ist das erforderlich? Ursache von Seung Bus Fehler?)
    914           nzb_local(ny+1:ny+gls,:) = nzb_local(0:gls-1,:)
     921!
     922!--       Exchange ghost points         
     923          CALL exchange_horiz_2d_int( nzb_local, nys, nyn, nxl, nxr, nbgp )
    915924       ENDIF
    916 
    917925!
    918926!--    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
    921929
    922930!
    … …  
    937945!--    nzb_s_outer:
    938946!--    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 + 1
    941           DO  i = 0, nx
     947       nzb_tmp = nzb_local
     948       DO  j = nys, nyn
     949          DO  i = nxl, nxr
    942950             nzb_tmp(j,i) = MAX( nzb_local(j,i-1), nzb_local(j,i),             &
    943951                                 nzb_local(j,i+1) )
    944952          ENDDO
    945953       ENDDO
     954       
     955       CALL exchange_horiz_2d_int( nzb_tmp, nys, nyn, nxl, nxr, nbgp )
     956       
    946957       DO  i = nxl, nxr
    947958          DO  j = nys, nyn
    … …  
    969980!--    nzb_u_inner:
    970981!--    extend nzb_local rightwards only
    971        nzb_tmp = nzb_local(-nbgp:ny+nbgp,-nbgp:nx+nbgp)
    972        DO  j = -1, ny + 1
    973           DO  i = 0, nx + 1
     982       nzb_tmp = nzb_local
     983       DO  j = nys, nyn
     984          DO  i = nxl, nxr
    974985             nzb_tmp(j,i) = MAX( nzb_local(j,i-1), nzb_local(j,i) )
    975986          ENDDO
    976987       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
    979992!
    980993!--    nzb_u_outer:
    … …  
    10011014!--    nzb_v_inner:
    10021015!--    extend nzb_local northwards only
    1003        nzb_tmp = nzb_local(-nbgp:ny+nbgp,-nbgp:nx+nbgp)
    1004        DO  i = -1, nx + 1
    1005           DO  j = 0, ny + 1
     1016       nzb_tmp = nzb_local
     1017       DO  i = nxl, nxr
     1018          DO  j = nys, nyn
    10061019             nzb_tmp(j,i) = MAX( nzb_local(j-1,i), nzb_local(j,i) )
    10071020          ENDDO
    10081021       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
    10101025
    10111026!
    … …  
    10351050!--    Since nzb_s_inner and nzb_w_inner are derived directly from nzb_local
    10361051!--    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 )
    10431058
    10441059!
    10451060!--    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))
    10551073          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
    10581106!
    10591107!--    Set flag arrays to be used for masking of grid points
    … …  
    10681116
    10691117    ENDIF
     1118!
     1119!-- Deallocate temporary array, as it might be reused for different
     1120!-- grid-levels further below.
     1121    DEALLOCATE( nzb_tmp )
    10701122
    10711123!
    … …  
    11901242       ENDDO
    11911243    ENDDO
    1192 
    11931244!
    11941245!-- Calculate wall flag arrays for the multigrid method.
    … …  
    11961247!-- version.
    11971248    IF ( psolver == 'multigrid_noopt' )  THEN
    1198 !
    1199 !--    Gridpoint increment of the current level
     1249
     1250!
     1251!--    Gridpoint increment of the current level.
    12001252       inc = 1
    1201 
    12021253       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
    12031257
    12041258          nxl_l = nxl_mg(l)
    … …  
    12071261          nyn_l = nyn_mg(l)
    12081262          nzt_l = nzt_mg(l)
    1209 
    12101263!
    12111264!--       Assign the flag level to be calculated
    … …  
    12511304          IF ( .NOT. masking_method )  THEN
    12521305
     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_l-2:nyn_l+2,nxl_l-2: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_l-1 / nyn_l+1 / nxl_l-1 / nxr_l+1. The alternative would be to
     1322!--          exchange 3D-INTEGER 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 non-cyclic 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                       
    12531347             DO  i = nxl_l-1, nxr_l+1
    12541348                DO  j = nys_l-1, nyn_l+1
    1255                    DO  k = nzb, nzt_l+1
    1256                          
     1349                   DO  k = nzb, nzt_l+1     
    12571350!
    12581351!--                   Inside/outside building (inside building does not need
    12591352!--                   further tests for walls)
    1260                       IF ( k*inc <= nzb_local(j*inc,i*inc) )  THEN
     1353                      IF ( k*inc <= nzb_tmp(j,i) )  THEN
    12611354
    12621355                         flags(k,j,i) = IBSET( flags(k,j,i), 6 )
    … …  
    12651358!
    12661359!--                      Bottom wall
    1267                          IF ( (k-1)*inc <= nzb_local(j*inc,i*inc) )  THEN
     1360                         IF ( (k-1)*inc <= nzb_tmp(j,i) )  THEN
    12681361                            flags(k,j,i) = IBSET( flags(k,j,i), 0 )
    12691362                         ENDIF
    12701363!
    12711364!--                      South wall
    1272                          IF ( k*inc <= nzb_local((j-1)*inc,i*inc) )  THEN
     1365                         IF ( k*inc <= nzb_tmp(j-1,i) )  THEN
    12731366                            flags(k,j,i) = IBSET( flags(k,j,i), 2 )
    12741367                         ENDIF
    12751368!
    12761369!--                      North wall
    1277                          IF ( k*inc <= nzb_local((j+1)*inc,i*inc) )  THEN
     1370                         IF ( k*inc <= nzb_tmp(j+1,i) )  THEN
    12781371                            flags(k,j,i) = IBSET( flags(k,j,i), 3 )
    12791372                         ENDIF
    12801373!
    12811374!--                      Left wall
    1282                          IF ( k*inc <= nzb_local(j*inc,(i-1)*inc) )  THEN
     1375                         IF ( k*inc <= nzb_tmp(j,i-1) )  THEN
    12831376                            flags(k,j,i) = IBSET( flags(k,j,i), 4 )
    12841377                         ENDIF
    12851378!
    12861379!--                      Right wall
    1287                          IF ( k*inc <= nzb_local(j*inc,(i+1)*inc) )  THEN
     1380                         IF ( k*inc <= nzb_tmp(j,i+1) )  THEN
    12881381                            flags(k,j,i) = IBSET( flags(k,j,i), 5 )
    12891382                         ENDIF
    … …  
    12951388             ENDDO
    12961389
     1390             DEALLOCATE( nzb_tmp )
     1391
    12971392          ENDIF
    12981393
    … …  
    13001395
    13011396       ENDDO
    1302 
     1397!
     1398!--    Reset grid_level to "normal" grid
     1399       grid_level = 0
     1400       
    13031401    ENDIF
    13041402!
    … …  
    14641562
    14651563    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 )
    14671565
    14681566
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |