Changeset 114 for palm/trunk/SOURCE


Ignore:
Timestamp:
Oct 10, 2007 12:03:15 AM (17 years ago)
Author:
raasch
Message:

preliminary updates for implementing buildings in poismg

Location:
palm/trunk/SOURCE
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/CURRENT_MODIFICATIONS

    r110 r114  
    11New:
    22---
     3
     4Pressure boundary conditions for vertical walls added to the multigrid solver.
     5They are applied using new wall flag arrays (wall_flags_..) which are defined
     6for each grid level. New argument gls added to routine user_init_grid
     7(user_interface).
     8
     9init_grid, init_pegrid, modules, user_interface
    310
    411
     
    916Errors:
    1017------
     18
     19Bugfix: pleft/pright changed to pnorth/psouth in sendrecv of particle tail
     20numbers along y (advec_particles)
     21
     22Bugfix: model_string needed a default value (combine_plot_fields)
     23
     24advec_particles, combine_plot_fields
  • palm/trunk/SOURCE/advec_particles.f90

    r110 r114  
    44! Actual revisions:
    55! -----------------
     6! Bugfix: pleft/pright changed to pnorth/psouth in sendrecv of particle tail
     7! numbers along y
    68! TEST: PRINT statements on unit 9 (commented out)
    79!
     
    27732775          IF ( use_particle_tails )  THEN
    27742776
    2775              CALL MPI_SENDRECV( trspt_count,      1, MPI_INTEGER, pleft, 0, &
    2776                                 trnpt_count_recv, 1, MPI_INTEGER, pright, 0, &
     2777             CALL MPI_SENDRECV( trspt_count,      1, MPI_INTEGER, psouth, 0, &
     2778                                trnpt_count_recv, 1, MPI_INTEGER, pnorth, 0, &
    27772779                                comm2d, status, ierr )
    27782780
     
    28502852          IF ( use_particle_tails )  THEN
    28512853
    2852              CALL MPI_SENDRECV( trnpt_count,      1, MPI_INTEGER, pright, 0, &
    2853                                 trspt_count_recv, 1, MPI_INTEGER, pleft, 0, &
     2854             CALL MPI_SENDRECV( trnpt_count,      1, MPI_INTEGER, pnorth, 0, &
     2855                                trspt_count_recv, 1, MPI_INTEGER, psouth, 0, &
    28542856                                comm2d, status, ierr )
    28552857
  • palm/trunk/SOURCE/check_parameters.f90

    r110 r114  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Multigrid solver allows topography
    77!
    88! Former revisions:
     
    286286          WRITE( action, '(A,A)' )  'timestep_scheme = ', timestep_scheme
    287287       ENDIF
    288        IF ( psolver == 'multigrid'  .OR.  psolver == 'sor' )  THEN
     288       IF ( psolver == 'sor' )  THEN
    289289          WRITE( action, '(A,A)' )  'psolver = ', psolver
    290290       ENDIF
  • palm/trunk/SOURCE/init_grid.f90

    r98 r114  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Calculation of wall flag arrays
    77!
    88! Former revisions:
     
    4444    IMPLICIT NONE
    4545
    46     INTEGER ::  bh, blx, bly, bxl, bxr, byn, bys, i, i_center, j, j_center, k, &
    47                 l, nzb_si, nzt_l, vi
     46    INTEGER ::  bh, blx, bly, bxl, bxr, byn, bys, gls, i, inc, i_center, j, &
     47                j_center, k, l, nxl_l, nxr_l, nyn_l, nys_l, nzb_si, nzt_l, vi
    4848
    4949    INTEGER, DIMENSION(:), ALLOCATABLE   ::  vertical_influence
     
    217217!
    218218!-- Allocate outer and inner index arrays for topography and set
    219 !-- defaults
    220     ALLOCATE( corner_nl(nys:nyn,nxl:nxr), corner_nr(nys:nyn,nxl:nxr), &
    221               corner_sl(nys:nyn,nxl:nxr), corner_sr(nys:nyn,nxl:nxr), &
    222               nzb_local(-1:ny+1,-1:nx+1), nzb_tmp(-1:ny+1,-1:nx+1),   &
    223               wall_l(nys:nyn,nxl:nxr), wall_n(nys:nyn,nxl:nxr),       &
     219!-- defaults.
     220!-- nzb_local has to contain additional layers of ghost points for calculating
     221!-- the flag arrays needed for the multigrid method
     222    gls = 2**( maximum_grid_level )
     223    ALLOCATE( corner_nl(nys:nyn,nxl:nxr), corner_nr(nys:nyn,nxl:nxr),       &
     224              corner_sl(nys:nyn,nxl:nxr), corner_sr(nys:nyn,nxl:nxr),       &
     225              nzb_local(-gls:ny+gls,-gls:nx+gls), nzb_tmp(-1:ny+1,-1:nx+1), &
     226              wall_l(nys:nyn,nxl:nxr), wall_n(nys:nyn,nxl:nxr),             &
    224227              wall_r(nys:nyn,nxl:nxr), wall_s(nys:nyn,nxl:nxr) )
    225228    ALLOCATE( fwxm(nys-1:nyn+1,nxl-1:nxr+1), fwxp(nys-1:nyn+1,nxl-1:nxr+1), &
     
    388391             ENDDO
    389392          ENDDO
    390           nzb_local(-1,0:nx)   = nzb_local(ny,0:nx)
    391           nzb_local(ny+1,0:nx) = nzb_local(0,0:nx)
    392           nzb_local(:,-1)      = nzb_local(:,nx)
    393           nzb_local(:,nx+1)    = nzb_local(:,0)
     393!
     394!--       Add cyclic boundaries (additional layers are for calculating flag
     395!--       arrays needed for the multigrid sover)
     396          nzb_local(-gls:-1,0:nx)     = nzb_local(ny-gls+1:ny,0:nx)
     397          nzb_local(ny+1:ny+gls,0:nx) = nzb_local(0:gls-1,0:nx)
     398          nzb_local(:,-gls:-1)        = nzb_local(:,nx-gls+1:nx)
     399          nzb_local(:,nx+1:nx+gls)    = nzb_local(:,0:gls-1)
    394400     
    395401          GOTO 12
     
    413419!--       case in the user interface. There, the subroutine user_init_grid
    414420!--       checks which of these two conditions applies.
    415           CALL user_init_grid( nzb_local )
     421          CALL user_init_grid( gls, nzb_local )
    416422
    417423    END SELECT
     424
     425!
     426!-- Test output of nzb_local -1:ny+1,-1:nx+1
     427    WRITE (9,*)  '*** nzb_local ***'
     428    DO  j = ny+1, -1, -1
     429       WRITE (9,'(194(1X,I2))')  ( nzb_local(j,i), i = -1, nx+1 )
     430    ENDDO
    418431
    419432!
     
    478491!--    nzb_s_outer:
    479492!--    extend nzb_local east-/westwards first, then north-/southwards
    480        nzb_tmp = nzb_local
     493       nzb_tmp = nzb_local(-1:ny+1,-1:nx+1)
    481494       DO  j = -1, ny + 1
    482495          DO  i = 0, nx
     
    510523!--    nzb_u_inner:
    511524!--    extend nzb_local rightwards only
    512        nzb_tmp = nzb_local
     525       nzb_tmp = nzb_local(-1:ny+1,-1:nx+1)
    513526       DO  j = -1, ny + 1
    514527          DO  i = 0, nx + 1
     
    542555!--    nzb_v_inner:
    543556!--    extend nzb_local northwards only
    544        nzb_tmp = nzb_local
     557       nzb_tmp = nzb_local(-1:ny+1,-1:nx+1)
    545558       DO  i = -1, nx + 1
    546559          DO  j = 0, ny + 1
     
    728741
    729742!
     743!-- Calculate wall flag arrays for the multigrid method
     744    IF ( psolver == 'multigrid' )  THEN
     745!
     746!--    Gridpoint increment of the current level
     747       inc = 1
     748
     749       DO  l = maximum_grid_level, 1 , -1
     750
     751          nxl_l = nxl_mg(l)
     752          nxr_l = nxr_mg(l)
     753          nys_l = nys_mg(l)
     754          nyn_l = nyn_mg(l)
     755          nzt_l = nzt_mg(l)
     756
     757!
     758!--       Assign the flag level to be calculated
     759          SELECT CASE ( l )
     760             CASE ( 1 )
     761                flags => wall_flags_1
     762             CASE ( 2 )
     763                flags => wall_flags_2
     764             CASE ( 3 )
     765                flags => wall_flags_3
     766             CASE ( 4 )
     767                flags => wall_flags_4
     768             CASE ( 5 )
     769                flags => wall_flags_5
     770             CASE ( 6 )
     771                flags => wall_flags_6
     772             CASE ( 7 )
     773                flags => wall_flags_7
     774             CASE ( 8 )
     775                flags => wall_flags_8
     776             CASE ( 9 )
     777                flags => wall_flags_9
     778             CASE ( 10 )
     779                flags => wall_flags_10
     780          END SELECT
     781
     782!
     783!--       Depending on the grid level, set the respective bits in case of
     784!--       neighbouring walls
     785!--       Bit 0:  wall to the bottom
     786!--       Bit 1:  wall to the top (not realized in remaining PALM code so far)
     787!--       Bit 2:  wall to the south
     788!--       Bit 3:  wall to the north
     789!--       Bit 4:  wall to the left
     790!--       Bit 5:  wall to the right
     791!--       Bit 6:  inside / outside building (1/0)
     792
     793          flags = 0
     794
     795          DO  i = nxl_l-1, nxr_l+1
     796             DO  j = nys_l-1, nyn_l+1
     797                DO  k = nzb, nzt_l+1
     798                         
     799!
     800!--                Inside/outside building (inside building does not need
     801!--                further tests for walls)
     802                   IF ( k*inc <= nzb_local(j*inc,i*inc) )  THEN
     803
     804                      flags(k,j,i) = IBSET( flags(k,j,i), 6 )
     805
     806                   ELSE
     807!
     808!--                   Bottom wall
     809                      IF ( (k-1)*inc <= nzb_local(j*inc,i*inc) )  THEN
     810                         flags(k,j,i) = IBSET( flags(k,j,i), 0 )
     811                      ENDIF
     812!
     813!--                   South wall
     814                      IF ( k*inc <= nzb_local((j-1)*inc,i*inc) )  THEN
     815                         flags(k,j,i) = IBSET( flags(k,j,i), 2 )
     816                      ENDIF
     817!
     818!--                   North wall
     819                      IF ( k*inc <= nzb_local((j+1)*inc,i*inc) )  THEN
     820                         flags(k,j,i) = IBSET( flags(k,j,i), 3 )
     821                      ENDIF
     822!
     823!--                   Left wall
     824                      IF ( k*inc <= nzb_local(j*inc,(i-1)*inc) )  THEN
     825                         flags(k,j,i) = IBSET( flags(k,j,i), 4 )
     826                      ENDIF
     827!
     828!--                   Right wall
     829                      IF ( k*inc <= nzb_local(j*inc,(i+1)*inc) )  THEN
     830                         flags(k,j,i) = IBSET( flags(k,j,i), 5 )
     831                      ENDIF
     832
     833                   ENDIF
     834                           
     835                ENDDO
     836             ENDDO
     837          ENDDO
     838
     839!
     840!--       Test output of flag arrays
     841          i = nxl_l
     842          WRITE (9,*)  ' '
     843          WRITE (9,*)  '*** mg level ', l, ' ***', mg_switch_to_pe0_level
     844          WRITE (9,*)  '    inc=', inc, '  i =', nxl_l
     845          WRITE (9,*)  '    nxl_l',nxl_l,' nxr_l=',nxr_l,' nys_l=',nys_l,' nyn_l=',nyn_l
     846          DO  k = nzt_l+1, nzb, -1
     847             WRITE (9,'(194(1X,I2))')  ( flags(k,j,i), j = nys_l-1, nyn_l+1 )
     848          ENDDO
     849
     850          inc = inc * 2
     851
     852       ENDDO
     853
     854    ENDIF
     855
     856!
    730857!-- In case of topography: limit near-wall mixing length l_wall further:
    731858!-- Go through all points of the subdomain one by one and look for the closest
  • palm/trunk/SOURCE/init_pegrid.f90

    r110 r114  
    44! Actual revisions:
    55! -----------------
     6! Allocation of wall flag arrays for multigrid solver
    67! TEST OUTPUT (TO BE REMOVED) logging mpi2 ierr values
    78!
     
    915916    ENDIF
    916917
     918!
     919!-- Allocate wall flag arrays used in the multigrid solver
     920    IF ( psolver == 'multigrid' )  THEN
     921
     922       DO  i = maximum_grid_level, 1, -1
     923
     924           SELECT CASE ( i )
     925
     926              CASE ( 1 )
     927                 ALLOCATE( wall_flags_1(nzb:nzt_mg(i)+1,         &
     928                                        nys_mg(i)-1:nyn_mg(i)+1, &
     929                                        nxl_mg(i)-1:nxr_mg(i)+1) )
     930
     931              CASE ( 2 )
     932                 ALLOCATE( wall_flags_2(nzb:nzt_mg(i)+1,         &
     933                                        nys_mg(i)-1:nyn_mg(i)+1, &
     934                                        nxl_mg(i)-1:nxr_mg(i)+1) )
     935
     936              CASE ( 3 )
     937                 ALLOCATE( wall_flags_3(nzb:nzt_mg(i)+1,         &
     938                                        nys_mg(i)-1:nyn_mg(i)+1, &
     939                                        nxl_mg(i)-1:nxr_mg(i)+1) )
     940
     941              CASE ( 4 )
     942                 ALLOCATE( wall_flags_4(nzb:nzt_mg(i)+1,         &
     943                                        nys_mg(i)-1:nyn_mg(i)+1, &
     944                                        nxl_mg(i)-1:nxr_mg(i)+1) )
     945
     946              CASE ( 5 )
     947                 ALLOCATE( wall_flags_5(nzb:nzt_mg(i)+1,         &
     948                                        nys_mg(i)-1:nyn_mg(i)+1, &
     949                                        nxl_mg(i)-1:nxr_mg(i)+1) )
     950
     951              CASE ( 6 )
     952                 ALLOCATE( wall_flags_6(nzb:nzt_mg(i)+1,         &
     953                                        nys_mg(i)-1:nyn_mg(i)+1, &
     954                                        nxl_mg(i)-1:nxr_mg(i)+1) )
     955
     956              CASE ( 7 )
     957                 ALLOCATE( wall_flags_7(nzb:nzt_mg(i)+1,         &
     958                                        nys_mg(i)-1:nyn_mg(i)+1, &
     959                                        nxl_mg(i)-1:nxr_mg(i)+1) )
     960
     961              CASE ( 8 )
     962                 ALLOCATE( wall_flags_8(nzb:nzt_mg(i)+1,         &
     963                                        nys_mg(i)-1:nyn_mg(i)+1, &
     964                                        nxl_mg(i)-1:nxr_mg(i)+1) )
     965
     966              CASE ( 9 )
     967                 ALLOCATE( wall_flags_9(nzb:nzt_mg(i)+1,         &
     968                                        nys_mg(i)-1:nyn_mg(i)+1, &
     969                                        nxl_mg(i)-1:nxr_mg(i)+1) )
     970
     971              CASE ( 10 )
     972                 ALLOCATE( wall_flags_10(nzb:nzt_mg(i)+1,        &
     973                                        nys_mg(i)-1:nyn_mg(i)+1, &
     974                                        nxl_mg(i)-1:nxr_mg(i)+1) )
     975
     976              CASE DEFAULT
     977                 IF ( myid == 0 )  PRINT*, '+++ init_pegrid: more than 10 ', &
     978                                           ' multigrid levels'
     979                 CALL local_stop
     980
     981          END SELECT
     982
     983       ENDDO
     984
     985    ENDIF
     986
    917987 END SUBROUTINE init_pegrid
  • palm/trunk/SOURCE/modules.f90

    r110 r114  
    55! Actual revisions:
    66! -----------------
    7 !
     7! +flags, wall_flags_1..10
    88!
    99! Former revisions:
     
    589589                nzb_v_outer, nzb_w_inner, nzb_w_outer, nzb_2d
    590590
     591    INTEGER, DIMENSION(:,:,:), POINTER ::  flags
     592
     593    INTEGER, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  wall_flags_1,           &
     594                wall_flags_2, wall_flags_3, wall_flags_4, wall_flags_5,        &
     595                wall_flags_6, wall_flags_7, wall_flags_8, wall_flags_9,        &
     596                wall_flags_10
     597
    591598    SAVE
    592599
  • palm/trunk/SOURCE/palm.f90

    r110 r114  
    6666    INTEGER           ::  i, run_description_header_i(80)
    6767
    68     version = 'PALM 3.4'
     68    version = 'PALM 3.4a'
    6969
    7070#if defined( __parallel )
  • palm/trunk/SOURCE/poismg.f90

    r77 r114  
    88! Actual revisions:
    99! -----------------
    10 !
     10! Boundary conditions at walls are implicitly set using flag arrays. Only
     11! Neumann BC is allowed. Upper walls are still not realized.
     12! Bottom and top BCs for array f_mg in restrict removed because boundary
     13! values are not needed (right hand side of SOR iteration)
    1114!
    1215! Former revisions:
     
    150153    l = grid_level
    151154
     155!
     156!-- Choose flag array of this level
     157    SELECT CASE ( l )
     158       CASE ( 1 )
     159          flags => wall_flags_1
     160       CASE ( 2 )
     161          flags => wall_flags_2
     162       CASE ( 3 )
     163          flags => wall_flags_3
     164       CASE ( 4 )
     165          flags => wall_flags_4
     166       CASE ( 5 )
     167          flags => wall_flags_5
     168       CASE ( 6 )
     169          flags => wall_flags_6
     170       CASE ( 7 )
     171          flags => wall_flags_7
     172       CASE ( 8 )
     173          flags => wall_flags_8
     174       CASE ( 9 )
     175          flags => wall_flags_9
     176       CASE ( 10 )
     177          flags => wall_flags_10
     178    END SELECT
     179
    152180!$OMP PARALLEL PRIVATE (i,j,k)
    153181!$OMP DO
     
    155183       DO  j = nys_mg(l), nyn_mg(l)
    156184          DO  k = nzb+1, nzt_mg(l)
    157              r(k,j,i) = f_mg(k,j,i)                                      &
    158                         - ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
    159                         - ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
    160                         - f2_mg(k,l) * p_mg(k+1,j,i)                     &
    161                         - f3_mg(k,l) * p_mg(k-1,j,i)                     &
     185             r(k,j,i) = f_mg(k,j,i)                                         &
     186                        - ddx2_mg(l) *                                      &
     187                            ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
     188                                          ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
     189                              p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
     190                                          ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
     191                        - ddy2_mg(l) *                                      &
     192                            ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
     193                                          ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
     194                              p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
     195                                          ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
     196                        - f2_mg(k,l) * p_mg(k+1,j,i)                        &
     197                        - f3_mg(k,l) *                                      &
     198                            ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
     199                                          ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
    162200                        + f1_mg(k,l) * p_mg(k,j,i)
     201!
     202!--          Residual within topography should be zero
     203             r(k,j,i) = r(k,j,i) * ( 1.0 - IBITS( flags(k,j,i), 6, 1 ) )
    163204          ENDDO
    164205       ENDDO
     
    181222
    182223!
    183 !-- Bottom and top boundary conditions
    184     IF ( ibc_p_b == 1 )  THEN
    185        r(nzb,:,: ) = r(nzb+1,:,:)
    186     ELSE
    187        r(nzb,:,: ) = 0.0
    188     ENDIF
    189 
     224!-- Top boundary condition
     225!-- A Neumann boundary condition for r is implicitly set in routine restrict
    190226    IF ( ibc_p_t == 1 )  THEN
    191227       r(nzt_mg(l)+1,:,: ) = r(nzt_mg(l),:,:)
     
    217253    INTEGER ::  i, ic, j, jc, k, kc, l
    218254
     255    REAL ::  rkjim, rkjip, rkjmi, rkjmim, rkjmip, rkjpi, rkjpim, rkjpip,       &
     256             rkmji, rkmjim, rkmjip, rkmjmi, rkmjmim, rkmjmip, rkmjpi, rkmjpim, &
     257             rkmjpip
     258
    219259    REAL, DIMENSION(nzb:nzt_mg(grid_level)+1,                            &
    220260                    nys_mg(grid_level)-1:nyn_mg(grid_level)+1,           &
     
    229269    l = grid_level
    230270
     271!
     272!-- Choose flag array of the upper level
     273    SELECT CASE ( l )
     274       CASE ( 1 )
     275          flags => wall_flags_1
     276       CASE ( 2 )
     277          flags => wall_flags_2
     278       CASE ( 3 )
     279          flags => wall_flags_3
     280       CASE ( 4 )
     281          flags => wall_flags_4
     282       CASE ( 5 )
     283          flags => wall_flags_5
     284       CASE ( 6 )
     285          flags => wall_flags_6
     286       CASE ( 7 )
     287          flags => wall_flags_7
     288       CASE ( 8 )
     289          flags => wall_flags_8
     290       CASE ( 9 )
     291          flags => wall_flags_9
     292       CASE ( 10 )
     293          flags => wall_flags_10
     294    END SELECT
     295   
    231296!$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc)
    232297!$OMP DO
     
    237302          DO  kc = nzb+1, nzt_mg(l)
    238303             k = 2*kc-1
     304!
     305!--          Use implicit Neumann BCs if the respective gridpoint is inside
     306!--          the building
     307             rkjim   = r(k,j,i-1)     + IBITS( flags(k,j,i-1), 6, 1 ) *     &
     308                                        ( r(k,j,i) - r(k,j,i-1) )
     309             rkjip   = r(k,j,i+1)     + IBITS( flags(k,j,i+1), 6, 1 ) *     &
     310                                        ( r(k,j,i) - r(k,j,i+1) )
     311             rkjpi   = r(k,j+1,i)     + IBITS( flags(k,j+1,i), 6, 1 ) *     &
     312                                        ( r(k,j,i) - r(k,j+1,i) )
     313             rkjmi   = r(k,j-1,i)     + IBITS( flags(k,j-1,i), 6, 1 ) *     &
     314                                        ( r(k,j,i) - r(k,j-1,i) )
     315             rkjmim  = r(k,j-1,i-1)   + IBITS( flags(k,j-1,i-1), 6, 1 ) *   &
     316                                        ( r(k,j,i) - r(k,j-1,i-1) )
     317             rkjpim  = r(k,j+1,i-1)   + IBITS( flags(k,j+1,i-1), 6, 1 ) *   &
     318                                        ( r(k,j,i) - r(k,j+1,i-1) )
     319             rkjmip  = r(k,j-1,i+1)   + IBITS( flags(k,j-1,i+1), 6, 1 ) *   &
     320                                        ( r(k,j,i) - r(k,j-1,i+1) )
     321             rkjpip  = r(k,j+1,i+1)   + IBITS( flags(k,j+1,i+1), 6, 1 ) *   &
     322                                        ( r(k,j,i) - r(k,j+1,i+1) )
     323             rkmji   = r(k-1,j,i)     + IBITS( flags(k-1,j,i), 6, 1 ) *     &
     324                                        ( r(k,j,i) - r(k-1,j,i) )
     325             rkmjim  = r(k-1,j,i-1)   + IBITS( flags(k-1,j,i-1), 6, 1 ) *   &
     326                                        ( r(k,j,i) - r(k-1,j,i-1) )
     327             rkmjip  = r(k-1,j,i+1)   + IBITS( flags(k-1,j,i+1), 6, 1 ) *   &
     328                                        ( r(k,j,i) - r(k-1,j,i+1) )
     329             rkmjpi  = r(k-1,j+1,i)   + IBITS( flags(k-1,j+1,i), 6, 1 ) *   &
     330                                        ( r(k,j,i) - r(k-1,j+1,i) )
     331             rkmjmi  = r(k-1,j-1,i)   + IBITS( flags(k-1,j-1,i), 6, 1 ) *   &
     332                                        ( r(k,j,i) - r(k-1,j-1,i) )
     333             rkmjmim = r(k-1,j-1,i-1) + IBITS( flags(k-1,j-1,i-1), 6, 1 ) * &
     334                                        ( r(k,j,i) - r(k-1,j-1,i-1) )
     335             rkmjpim = r(k-1,j+1,i-1) + IBITS( flags(k-1,j+1,i-1), 6, 1 ) * &
     336                                        ( r(k,j,i) - r(k-1,j+1,i-1) )
     337             rkmjmip = r(k-1,j-1,i+1) + IBITS( flags(k-1,j-1,i+1), 6, 1 ) * &
     338                                        ( r(k,j,i) - r(k-1,j-1,i+1) )
     339             rkmjpip = r(k-1,j+1,i+1) + IBITS( flags(k-1,j+1,i+1), 6, 1 ) * &
     340                                        ( r(k,j,i) - r(k-1,j+1,i+1) )
     341
    239342             f_mg(kc,jc,ic) = 1.0 / 64.0 * (                            &
    240343                              8.0 * r(k,j,i)                            &
    241                             + 4.0 * ( r(k,j,i-1)     + r(k,j,i+1)     + &
    242                                       r(k,j+1,i)     + r(k,j-1,i)     ) &
    243                             + 2.0 * ( r(k,j-1,i-1)   + r(k,j+1,i-1)   + &
    244                                       r(k,j-1,i+1)   + r(k,j+1,i+1)   ) &
    245                             + 4.0 * r(k-1,j,i)                          &
    246                             + 2.0 * ( r(k-1,j,i-1)   + r(k-1,j,i+1)   + &
    247                                       r(k-1,j+1,i)   + r(k-1,j-1,i)   ) &
    248                             +       ( r(k-1,j-1,i-1) + r(k-1,j+1,i-1) + &
    249                                       r(k-1,j-1,i+1) + r(k-1,j+1,i+1) ) &
     344                            + 4.0 * ( rkjim   + rkjip   +              &
     345                                      rkjpi   + rkjmi   )              &
     346                            + 2.0 * ( rkjmim  + rkjpim  +              &
     347                                      rkjmip  + rkjpip  )              &
     348                            + 4.0 * rkmji                               &
     349                            + 2.0 * ( rkmjim  + rkmjim  +              &
     350                                      rkmjpi  + rkmjmi  )              &
     351                            +       ( rkmjmim + rkmjpim +              &
     352                                      rkmjmip + rkmjpip )              &
    250353                            + 4.0 * r(k+1,j,i)                          &
    251354                            + 2.0 * ( r(k+1,j,i-1)   + r(k+1,j,i+1)   + &
     
    254357                                      r(k+1,j-1,i+1) + r(k+1,j+1,i+1) ) &
    255358                                           )
     359
     360!             f_mg(kc,jc,ic) = 1.0 / 64.0 * (                            &
     361!                              8.0 * r(k,j,i)                            &
     362!                            + 4.0 * ( r(k,j,i-1)     + r(k,j,i+1)     + &
     363!                                      r(k,j+1,i)     + r(k,j-1,i)     ) &
     364!                            + 2.0 * ( r(k,j-1,i-1)   + r(k,j+1,i-1)   + &
     365!                                      r(k,j-1,i+1)   + r(k,j+1,i+1)   ) &
     366!                            + 4.0 * r(k-1,j,i)                          &
     367!                            + 2.0 * ( r(k-1,j,i-1)   + r(k-1,j,i+1)   + &
     368!                                      r(k-1,j+1,i)   + r(k-1,j-1,i)   ) &
     369!                            +       ( r(k-1,j-1,i-1) + r(k-1,j+1,i-1) + &
     370!                                      r(k-1,j-1,i+1) + r(k-1,j+1,i+1) ) &
     371!                            + 4.0 * r(k+1,j,i)                          &
     372!                            + 2.0 * ( r(k+1,j,i-1)   + r(k+1,j,i+1)   + &
     373!                                      r(k+1,j+1,i)   + r(k+1,j-1,i)   ) &
     374!                            +       ( r(k+1,j-1,i-1) + r(k+1,j+1,i-1) + &
     375!                                      r(k+1,j-1,i+1) + r(k+1,j+1,i+1) ) &
     376!                                           )
    256377          ENDDO
    257378       ENDDO
     
    275396!
    276397!-- Bottom and top boundary conditions
    277     IF ( ibc_p_b == 1 )  THEN
    278        f_mg(nzb,:,: ) = f_mg(nzb+1,:,:)
    279     ELSE
    280        f_mg(nzb,:,: ) = 0.0
    281     ENDIF
    282 
    283     IF ( ibc_p_t == 1 )  THEN
    284        f_mg(nzt_mg(l)+1,:,: ) = f_mg(nzt_mg(l),:,:)
    285     ELSE
    286        f_mg(nzt_mg(l)+1,:,: ) = 0.0
    287     ENDIF
     398!    IF ( ibc_p_b == 1 )  THEN
     399!       f_mg(nzb,:,: ) = f_mg(nzb+1,:,:)
     400!    ELSE
     401!       f_mg(nzb,:,: ) = 0.0
     402!    ENDIF
     403!
     404!    IF ( ibc_p_t == 1 )  THEN
     405!       f_mg(nzt_mg(l)+1,:,: ) = f_mg(nzt_mg(l),:,:)
     406!    ELSE
     407!       f_mg(nzt_mg(l)+1,:,: ) = 0.0
     408!    ENDIF
    288409
    289410
     
    412533    LOGICAL :: unroll
    413534
     535    REAL ::  wall_left, wall_north, wall_right, wall_south, wall_total, wall_top
     536
    414537    REAL, DIMENSION(nzb:nzt_mg(grid_level)+1,                                 &
    415538                    nys_mg(grid_level)-1:nyn_mg(grid_level)+1,                &
     
    418541
    419542    l = grid_level
     543
     544!
     545!-- Choose flag array of this level
     546    SELECT CASE ( l )
     547       CASE ( 1 )
     548          flags => wall_flags_1
     549       CASE ( 2 )
     550          flags => wall_flags_2
     551       CASE ( 3 )
     552          flags => wall_flags_3
     553       CASE ( 4 )
     554          flags => wall_flags_4
     555       CASE ( 5 )
     556          flags => wall_flags_5
     557       CASE ( 6 )
     558          flags => wall_flags_6
     559       CASE ( 7 )
     560          flags => wall_flags_7
     561       CASE ( 8 )
     562          flags => wall_flags_8
     563       CASE ( 9 )
     564          flags => wall_flags_9
     565       CASE ( 10 )
     566          flags => wall_flags_10
     567    END SELECT
    420568
    421569    unroll = ( MOD( nyn_mg(l)-nys_mg(l)+1, 4 ) == 0  .AND. &
     
    434582                DO  j = nys_mg(l) + 2 - colour, nyn_mg(l), 2
    435583                   DO  k = nzb+1, nzt_mg(l), 2
     584!                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
     585!                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
     586!                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
     587!                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
     588!                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
     589!                                                       )
     590
    436591                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
    437                                 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
    438                               + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
    439                               + f2_mg(k,l) * p_mg(k+1,j,i)                     &
    440                               + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
    441                                                        )
     592                             ddx2_mg(l) *                                      &
     593                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
     594                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
     595                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
     596                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
     597                           + ddy2_mg(l) *                                      &
     598                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
     599                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
     600                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
     601                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
     602                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
     603                           + f3_mg(k,l) *                                      &
     604                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
     605                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
     606                           - f_mg(k,j,i)               )
    442607                   ENDDO
    443608                ENDDO
     
    448613                   DO  k = nzb+1, nzt_mg(l), 2
    449614                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
    450                                 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
    451                               + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
    452                               + f2_mg(k,l) * p_mg(k+1,j,i)                     &
    453                               + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
    454                                                        )
     615                             ddx2_mg(l) *                                      &
     616                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
     617                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
     618                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
     619                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
     620                           + ddy2_mg(l) *                                      &
     621                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
     622                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
     623                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
     624                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
     625                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
     626                           + f3_mg(k,l) *                                      &
     627                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
     628                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
     629                           - f_mg(k,j,i)               )
    455630                   ENDDO
    456631                ENDDO
     
    461636                   DO  k = nzb+2, nzt_mg(l), 2
    462637                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
    463                                 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
    464                               + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
    465                               + f2_mg(k,l) * p_mg(k+1,j,i)                     &
    466                               + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
    467                                                        )
     638                             ddx2_mg(l) *                                      &
     639                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
     640                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
     641                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
     642                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
     643                           + ddy2_mg(l) *                                      &
     644                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
     645                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
     646                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
     647                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
     648                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
     649                           + f3_mg(k,l) *                                      &
     650                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
     651                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
     652                           - f_mg(k,j,i)               )
    468653                   ENDDO
    469654                ENDDO
     
    474659                   DO  k = nzb+2, nzt_mg(l), 2
    475660                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
    476                                 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
    477                               + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
    478                               + f2_mg(k,l) * p_mg(k+1,j,i)                     &
    479                               + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
    480                                                        )
     661                             ddx2_mg(l) *                                      &
     662                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
     663                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
     664                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
     665                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
     666                           + ddy2_mg(l) *                                      &
     667                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
     668                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
     669                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
     670                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
     671                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
     672                           + f3_mg(k,l) *                                      &
     673                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
     674                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
     675                           - f_mg(k,j,i)               )
    481676                   ENDDO
    482677                ENDDO
     
    496691                      j = jj
    497692                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
    498                                 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
    499                               + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
    500                               + f2_mg(k,l) * p_mg(k+1,j,i)                     &
    501                               + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
    502                                                        )
     693                             ddx2_mg(l) *                                      &
     694                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
     695                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
     696                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
     697                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
     698                           + ddy2_mg(l) *                                      &
     699                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
     700                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
     701                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
     702                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
     703                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
     704                           + f3_mg(k,l) *                                      &
     705                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
     706                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
     707                           - f_mg(k,j,i)               )
    503708                      j = jj+2
    504709                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
    505                                 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
    506                               + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
    507                               + f2_mg(k,l) * p_mg(k+1,j,i)                     &
    508                               + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
    509                                                        )
    510 !                      j = jj+4
    511 !                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
    512 !                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
    513 !                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
    514 !                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
    515 !                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
    516 !                                                    )
    517 !                      j = jj+6
    518 !                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
    519 !                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
    520 !                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
    521 !                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
    522 !                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
    523 !                                                       )
     710                             ddx2_mg(l) *                                      &
     711                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
     712                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
     713                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
     714                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
     715                           + ddy2_mg(l) *                                      &
     716                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
     717                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
     718                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
     719                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
     720                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
     721                           + f3_mg(k,l) *                                      &
     722                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
     723                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
     724                           - f_mg(k,j,i)               )
    524725                   ENDDO
    525726   
     
    529730                      j =jj
    530731                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
    531                                 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
    532                               + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
    533                               + f2_mg(k,l) * p_mg(k+1,j,i)                     &
    534                               + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
    535                                                        )
     732                             ddx2_mg(l) *                                      &
     733                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
     734                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
     735                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
     736                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
     737                           + ddy2_mg(l) *                                      &
     738                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
     739                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
     740                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
     741                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
     742                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
     743                           + f3_mg(k,l) *                                      &
     744                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
     745                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
     746                           - f_mg(k,j,i)               )
    536747                      j = jj+2
    537748                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
    538                                 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
    539                               + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
    540                               + f2_mg(k,l) * p_mg(k+1,j,i)                     &
    541                               + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
    542                                                        )
    543 !                      j = jj+4
    544 !                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
    545 !                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
    546 !                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
    547 !                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
    548 !                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
    549 !                                                       )
    550 !                      j = jj+6
    551 !                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
    552 !                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
    553 !                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
    554 !                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
    555 !                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
    556 !                                                       )
     749                             ddx2_mg(l) *                                      &
     750                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
     751                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
     752                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
     753                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
     754                           + ddy2_mg(l) *                                      &
     755                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
     756                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
     757                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
     758                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
     759                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
     760                           + f3_mg(k,l) *                                      &
     761                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
     762                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
     763                           - f_mg(k,j,i)               )
    557764                   ENDDO
    558765
     
    562769                      j =jj
    563770                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
    564                                 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
    565                               + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
    566                               + f2_mg(k,l) * p_mg(k+1,j,i)                     &
    567                               + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
    568                                                        )
     771                             ddx2_mg(l) *                                      &
     772                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
     773                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
     774                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
     775                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
     776                           + ddy2_mg(l) *                                      &
     777                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
     778                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
     779                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
     780                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
     781                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
     782                           + f3_mg(k,l) *                                      &
     783                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
     784                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
     785                           - f_mg(k,j,i)               )
    569786                      j = jj+2
    570787                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
    571                                 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
    572                               + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
    573                               + f2_mg(k,l) * p_mg(k+1,j,i)                     &
    574                               + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
    575                                                        )
    576 !                      j = jj+4
    577 !                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
    578 !                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
    579 !                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
    580 !                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
    581 !                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
    582 !                                                       )
    583 !                      j = jj+6
    584 !                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
    585 !                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
    586 !                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
    587 !                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
    588 !                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
    589 !                                                       )
     788                             ddx2_mg(l) *                                      &
     789                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
     790                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
     791                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
     792                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
     793                           + ddy2_mg(l) *                                      &
     794                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
     795                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
     796                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
     797                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
     798                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
     799                           + f3_mg(k,l) *                                      &
     800                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
     801                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
     802                           - f_mg(k,j,i)               )
    590803                   ENDDO
    591804
     
    595808                      j =jj
    596809                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
    597                                 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
    598                               + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
    599                               + f2_mg(k,l) * p_mg(k+1,j,i)                     &
    600                               + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
    601                                                        )
     810                             ddx2_mg(l) *                                      &
     811                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
     812                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
     813                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
     814                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
     815                           + ddy2_mg(l) *                                      &
     816                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
     817                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
     818                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
     819                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
     820                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
     821                           + f3_mg(k,l) *                                      &
     822                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
     823                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
     824                           - f_mg(k,j,i)               )
    602825                      j = jj+2
    603826                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
    604                                 ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
    605                               + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
    606                               + f2_mg(k,l) * p_mg(k+1,j,i)                     &
    607                               + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
    608                                                        )
    609 !                      j = jj+4
    610 !                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
    611 !                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
    612 !                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
    613 !                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
    614 !                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
    615 !                                                       )
    616 !                      j = jj+6
    617 !                      p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                       &
    618 !                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
    619 !                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
    620 !                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
    621 !                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
    622 !                                                       )
     827                             ddx2_mg(l) *                                      &
     828                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
     829                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
     830                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
     831                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
     832                           + ddy2_mg(l) *                                      &
     833                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
     834                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
     835                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
     836                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
     837                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
     838                           + f3_mg(k,l) *                                      &
     839                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
     840                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
     841                           - f_mg(k,j,i)               )
    623842                   ENDDO
    624843
     
    669888    ENDDO
    670889
     890!
     891!-- Set pressure within topography and at the topography surfaces
     892!$OMP PARALLEL PRIVATE (i,j,k,wall_left,wall_north,wall_right,wall_south,wall_top,wall_total)
     893!$OMP DO
     894    DO  i = nxl_mg(l), nxr_mg(l)
     895       DO  j = nys_mg(l), nyn_mg(l)
     896          DO  k = nzb, nzt_mg(l)
     897!
     898!--          First, set pressure inside topography to zero
     899             p_mg(k,j,i) = p_mg(k,j,i) * ( 1.0 - IBITS( flags(k,j,i), 6, 1 ) )
     900!
     901!--          Second, determine if the gridpoint inside topography is adjacent
     902!--          to a wall and set its value to a value given by the average of
     903!--          those values obtained from Neumann boundary condition
     904             wall_left  = IBITS( flags(k,j,i-1), 5, 1 )
     905             wall_right = IBITS( flags(k,j,i+1), 4, 1 )
     906             wall_south = IBITS( flags(k,j-1,i), 3, 1 )
     907             wall_north = IBITS( flags(k,j+1,i), 2, 1 )
     908             wall_top   = IBITS( flags(k+1,j,i), 0, 1 )
     909             wall_total = wall_left + wall_right + wall_south + wall_north + &
     910                          wall_top
     911
     912             IF ( wall_total > 0.0 )  THEN
     913                p_mg(k,j,i) = 1.0 / wall_total *                 &
     914                                  ( wall_left  * p_mg(k,j,i-1) + &
     915                                    wall_right * p_mg(k,j,i+1) + &
     916                                    wall_south * p_mg(k,j-1,i) + &
     917                                    wall_north * p_mg(k,j+1,i) + &
     918                                    wall_top   * p_mg(k+1,j,i) )
     919             ENDIF
     920          ENDDO
     921       ENDDO
     922    ENDDO
     923!$OMP END PARALLEL
     924
     925!
     926!-- One more time horizontal boundary conditions
     927    CALL exchange_horiz( p_mg )
    671928
    672929 END SUBROUTINE redblack
  • palm/trunk/SOURCE/user_interface.f90

    r110 r114  
    44! Actual revisions:
    55! -----------------
    6 !
     6! +argument gls in user_init_grid
    77!
    88! Former revisions:
     
    238238
    239239
    240  SUBROUTINE user_init_grid( nzb_local )
     240 SUBROUTINE user_init_grid( gls, nzb_local )
    241241
    242242!------------------------------------------------------------------------------!
     
    245245! ------------
    246246! Execution of user-defined grid initializing actions
     247! First argument gls contains the number of ghost layers, which is > 1 if the
     248! multigrid method for the pressure solver is used
    247249!------------------------------------------------------------------------------!
    248250
     
    253255    IMPLICIT NONE
    254256
    255     INTEGER, DIMENSION(-1:ny+1,-1:nx+1) ::  nzb_local
     257    INTEGER ::  gls
     258
     259    INTEGER, DIMENSION(-gls:ny+gls,-gls:nx+gls) ::  nzb_local
    256260
    257261!
Note: See TracChangeset for help on using the changeset viewer.