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

preliminary updates for implementing buildings in poismg

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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
Note: See TracChangeset for help on using the changeset viewer.