Changeset 114 for palm/trunk/SOURCE/init_grid.f90

Timestamp:

Oct 10, 2007 12:03:15 AM (17 years ago)

Author:

raasch

Message:

preliminary updates for implementing buildings in poismg

File:

• palm/trunk/SOURCE/init_grid.f90 (modified) (9 diffs)

palm/trunk/SOURCE/init_grid.f90

@@ -4 +4 @@
 ! Actual revisions:
 ! -----------------
-! 
+! Calculation of wall flag arrays
 !
 ! Former revisions:
@@ -44 +44 @@
     IMPLICIT NONE
 
-    INTEGER ::  bh, blx, bly, bxl, bxr, byn, bys, i, i_center, j, j_center, k, &
-                l, nzb_si, nzt_l, vi
+    INTEGER ::  bh, blx, bly, bxl, bxr, byn, bys, gls, i, inc, i_center, j, &
+                j_center, k, l, nxl_l, nxr_l, nyn_l, nys_l, nzb_si, nzt_l, vi
 
     INTEGER, DIMENSION(:), ALLOCATABLE   ::  vertical_influence
@@ -217 +217 @@
 !
 !-- Allocate outer and inner index arrays for topography and set
-!-- defaults
-    ALLOCATE( corner_nl(nys:nyn,nxl:nxr), corner_nr(nys:nyn,nxl:nxr), &
-              corner_sl(nys:nyn,nxl:nxr), corner_sr(nys:nyn,nxl:nxr), &
-              nzb_local(-1:ny+1,-1:nx+1), nzb_tmp(-1:ny+1,-1:nx+1),   &
-              wall_l(nys:nyn,nxl:nxr), wall_n(nys:nyn,nxl:nxr),       &
+!-- defaults.
+!-- nzb_local has to contain additional layers of ghost points for calculating
+!-- the flag arrays needed for the multigrid method
+    gls = 2**( maximum_grid_level )
+    ALLOCATE( corner_nl(nys:nyn,nxl:nxr), corner_nr(nys:nyn,nxl:nxr),       &
+              corner_sl(nys:nyn,nxl:nxr), corner_sr(nys:nyn,nxl:nxr),       &
+              nzb_local(-gls:ny+gls,-gls:nx+gls), nzb_tmp(-1:ny+1,-1:nx+1), &
+              wall_l(nys:nyn,nxl:nxr), wall_n(nys:nyn,nxl:nxr),             &
               wall_r(nys:nyn,nxl:nxr), wall_s(nys:nyn,nxl:nxr) )
     ALLOCATE( fwxm(nys-1:nyn+1,nxl-1:nxr+1), fwxp(nys-1:nyn+1,nxl-1:nxr+1), &
@@ -388 +391 @@
              ENDDO
           ENDDO
-          nzb_local(-1,0:nx)   = nzb_local(ny,0:nx)
-          nzb_local(ny+1,0:nx) = nzb_local(0,0:nx)
-          nzb_local(:,-1)      = nzb_local(:,nx)
-          nzb_local(:,nx+1)    = nzb_local(:,0)
+!
+!--       Add cyclic boundaries (additional layers are for calculating flag
+!--       arrays needed for the multigrid sover)
+          nzb_local(-gls:-1,0:nx)     = nzb_local(ny-gls+1:ny,0:nx)
+          nzb_local(ny+1:ny+gls,0:nx) = nzb_local(0:gls-1,0:nx)
+          nzb_local(:,-gls:-1)        = nzb_local(:,nx-gls+1:nx)
+          nzb_local(:,nx+1:nx+gls)    = nzb_local(:,0:gls-1)
       
           GOTO 12
@@ -413 +419 @@
 !--       case in the user interface. There, the subroutine user_init_grid 
 !--       checks which of these two conditions applies.
-          CALL user_init_grid( nzb_local )
+          CALL user_init_grid( gls, nzb_local )
 
     END SELECT
+
+!
+!-- Test output of nzb_local -1:ny+1,-1:nx+1
+    WRITE (9,*)  '*** nzb_local ***'
+    DO  j = ny+1, -1, -1
+       WRITE (9,'(194(1X,I2))')  ( nzb_local(j,i), i = -1, nx+1 )
+    ENDDO
 
 !
@@ -478 +491 @@
 !--    nzb_s_outer: 
 !--    extend nzb_local east-/westwards first, then north-/southwards
-       nzb_tmp = nzb_local
+       nzb_tmp = nzb_local(-1:ny+1,-1:nx+1)
        DO  j = -1, ny + 1
           DO  i = 0, nx
@@ -510 +523 @@
 !--    nzb_u_inner: 
 !--    extend nzb_local rightwards only
-       nzb_tmp = nzb_local
+       nzb_tmp = nzb_local(-1:ny+1,-1:nx+1)
        DO  j = -1, ny + 1
           DO  i = 0, nx + 1
@@ -542 +555 @@
 !--    nzb_v_inner:
 !--    extend nzb_local northwards only
-       nzb_tmp = nzb_local
+       nzb_tmp = nzb_local(-1:ny+1,-1:nx+1)
        DO  i = -1, nx + 1
           DO  j = 0, ny + 1
@@ -728 +741 @@
 
 !
+!-- Calculate wall flag arrays for the multigrid method
+    IF ( psolver == 'multigrid' )  THEN
+!
+!--    Gridpoint increment of the current level
+       inc = 1
+
+       DO  l = maximum_grid_level, 1 , -1
+
+          nxl_l = nxl_mg(l)
+          nxr_l = nxr_mg(l)
+          nys_l = nys_mg(l)
+          nyn_l = nyn_mg(l)
+          nzt_l = nzt_mg(l)
+
+!
+!--       Assign the flag level to be calculated
+          SELECT CASE ( l )
+             CASE ( 1 )
+                flags => wall_flags_1
+             CASE ( 2 )
+                flags => wall_flags_2
+             CASE ( 3 )
+                flags => wall_flags_3
+             CASE ( 4 )
+                flags => wall_flags_4
+             CASE ( 5 )
+                flags => wall_flags_5
+             CASE ( 6 )
+                flags => wall_flags_6
+             CASE ( 7 )
+                flags => wall_flags_7
+             CASE ( 8 )
+                flags => wall_flags_8
+             CASE ( 9 )
+                flags => wall_flags_9
+             CASE ( 10 )
+                flags => wall_flags_10
+          END SELECT
+
+!
+!--       Depending on the grid level, set the respective bits in case of
+!--       neighbouring walls
+!--       Bit 0:  wall to the bottom
+!--       Bit 1:  wall to the top (not realized in remaining PALM code so far)
+!--       Bit 2:  wall to the south
+!--       Bit 3:  wall to the north
+!--       Bit 4:  wall to the left
+!--       Bit 5:  wall to the right
+!--       Bit 6:  inside / outside building (1/0)
+
+          flags = 0
+
+          DO  i = nxl_l-1, nxr_l+1
+             DO  j = nys_l-1, nyn_l+1
+                DO  k = nzb, nzt_l+1
+                         
+!
+!--                Inside/outside building (inside building does not need
+!--                further tests for walls)
+                   IF ( k*inc <= nzb_local(j*inc,i*inc) )  THEN
+
+                      flags(k,j,i) = IBSET( flags(k,j,i), 6 )
+
+                   ELSE
+!
+!--                   Bottom wall
+                      IF ( (k-1)*inc <= nzb_local(j*inc,i*inc) )  THEN
+                         flags(k,j,i) = IBSET( flags(k,j,i), 0 )
+                      ENDIF
+!
+!--                   South wall
+                      IF ( k*inc <= nzb_local((j-1)*inc,i*inc) )  THEN
+                         flags(k,j,i) = IBSET( flags(k,j,i), 2 )
+                      ENDIF
+!
+!--                   North wall
+                      IF ( k*inc <= nzb_local((j+1)*inc,i*inc) )  THEN
+                         flags(k,j,i) = IBSET( flags(k,j,i), 3 )
+                      ENDIF
+!
+!--                   Left wall
+                      IF ( k*inc <= nzb_local(j*inc,(i-1)*inc) )  THEN
+                         flags(k,j,i) = IBSET( flags(k,j,i), 4 )
+                      ENDIF
+!
+!--                   Right wall
+                      IF ( k*inc <= nzb_local(j*inc,(i+1)*inc) )  THEN
+                         flags(k,j,i) = IBSET( flags(k,j,i), 5 )
+                      ENDIF
+
+                   ENDIF
+                            
+                ENDDO
+             ENDDO
+          ENDDO 
+
+!
+!--       Test output of flag arrays
+          i = nxl_l
+          WRITE (9,*)  ' '
+          WRITE (9,*)  '*** mg level ', l, ' ***', mg_switch_to_pe0_level
+          WRITE (9,*)  '    inc=', inc, '  i =', nxl_l
+          WRITE (9,*)  '    nxl_l',nxl_l,' nxr_l=',nxr_l,' nys_l=',nys_l,' nyn_l=',nyn_l
+          DO  k = nzt_l+1, nzb, -1
+             WRITE (9,'(194(1X,I2))')  ( flags(k,j,i), j = nys_l-1, nyn_l+1 )
+          ENDDO
+
+          inc = inc * 2
+
+       ENDDO
+
+    ENDIF
+
+!
 !-- In case of topography: limit near-wall mixing length l_wall further:
 !-- Go through all points of the subdomain one by one and look for the closest