Changeset 114 for palm/trunk


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
Files:
12 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/DOC/app/chapter_4.1.html

    r111 r114  
    1616 
    1717  <title>PALM chapter 4.1</title></head>
    18 
    1918<body>
    2019
     
    92609259
    92619260In case of ocean runs (see <a href="chapter_4.1.html#ocean">ocean</a>),
    9262 this parameter gives the velocity value at the sea surface, which is
     9261this parameter gives the geostrophic velocity value (i.e. the pressure gradient) at the sea surface, which is
    92639262at k=nzt. The profile is then constructed from the surface down to the
    92649263bottom of the model.<br>
     
    1016810167
    1016910168In case of ocean runs (see <a href="chapter_4.1.html#ocean">ocean</a>),
    10170 this parameter gives the velocity value at the sea surface, which is
     10169this parameter gives the geostrophic velocity value (i.e. the pressure gradient) at the sea surface, which is
    1017110170at k=nzt. The profile is then constructed from the surface down to the
    1017210171bottom of the model.</td>
  • palm/trunk/SCRIPTS/mrun

    r109 r114  
    200200 read_from_config=""
    201201 restart_run=false
    202  return_addres=$(nslookup `hostname` 2>&1 | grep "Address:" | tail -1 | awk '{print $2}')
     202# return_addres=$(nslookup `hostname` 2>&1 | grep "Address:" | tail -1 | awk '{print $2}')
    203203 if [[ $return_addres = 130.75.105.158 ]]
    204204 then
  • 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!
  • palm/trunk/UTIL/combine_plot_fields.f90

    • Property svn:keywords set to Id
    r108 r114  
    44! Actual revisions:
    55! -----------------
    6 ! loop for processing of output by coupled runs, id_string does not contain
    7 ! modus any longer
     6! Bugfix: model_string needed a default value
    87!
    98! Former revisions:
    109! -----------------
     10! $Id$
     11!
     12! Aug 07    Loop for processing of output by coupled runs, id_string does not
     13!           contain modus any longer
     14!
    1115! 18/01/06  Output of time-averaged data
    1216!
     
    106110
    107111    PRINT*, ''
     112    PRINT*, ''
    108113    PRINT*, '*** combine_plot_fields ***'
     114
     115!
     116!-- Find out if a coupled run has been carried out
    109117    INQUIRE( FILE='COUPLING_PORT_OPENED', EXIST=found )
    110118    IF ( found )  THEN
     
    115123       PRINT*, '    uncoupled run'
    116124    ENDIF
     125
     126!
     127!-- Do everything for each model
    117128    DO model = 1, models
     129!
     130!--    Set the model string used to identify the filenames
     131       model_string = ''
    118132       IF ( models == 2 )  THEN
    119           PRINT*, ''
    120133          PRINT*, ''
    121134          PRINT*, '*** combine_plot_fields ***'
    122135          IF ( model == 2 )  THEN
    123              model_string = "_O"
     136             model_string = '_O'
    124137             PRINT*, '    now combining ocean data'
    125138             PRINT*, '    ========================'
    126139          ELSE
    127              model_string = ""
    128140             PRINT*, '    now combining atmosphere data'
    129141             PRINT*, '    ============================='
Note: See TracChangeset for help on using the changeset viewer.