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