Changeset 1628


Ignore:
Timestamp:
Aug 26, 2015 4:54:34 PM (6 years ago)
Author:
suehring
Message:

Bugfix concerning wall_flags at left and south PE boundaries

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/advec_ws.f90

    r1586 r1628  
    493493          DO  k = nzb+1, nzb_max
    494494
    495              ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
    496              ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
    497              ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
     495             ibit5 = IBITS(wall_flags_0(k,j-1,i),5,1)
     496             ibit4 = IBITS(wall_flags_0(k,j-1,i),4,1)
     497             ibit3 = IBITS(wall_flags_0(k,j-1,i),3,1)
    498498
    499499             v_comp                  = v(k,j,i) - v_gtrans
     
    554554          DO  k = nzb+1, nzb_max
    555555
    556              ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
    557              ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
    558              ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
     556             ibit2 = IBITS(wall_flags_0(k,j,i-1),2,1)
     557             ibit1 = IBITS(wall_flags_0(k,j,i-1),1,1)
     558             ibit0 = IBITS(wall_flags_0(k,j,i-1),0,1)
    559559
    560560             u_comp                     = u(k,j,i) - u_gtrans
     
    874874!--       correction is needed to overcome numerical instabilities caused
    875875!--       by a not sufficient reduction of divergences near topography.
    876           div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx                   &
    877                         + ( v(k,j+1,i) - v(k,j,i)   ) * ddy                   &
    878                         + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
     876          div         =   ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 )             &
     877                          - u(k,j,i)   * ( IBITS(wall_flags_0(k,j,i-1),0,1)    &
     878                                         + IBITS(wall_flags_0(k,j,i-1),1,1)    &
     879                                         + IBITS(wall_flags_0(k,j,i-1),2,1)    &
     880                                         )                                     &
     881                          ) * ddx                                              &
     882                        + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )             &
     883                          - v(k,j,i)   * ( IBITS(wall_flags_0(k,j-1,i),3,1)    &
     884                                         + IBITS(wall_flags_0(k,j-1,i),4,1)    &
     885                                         + IBITS(wall_flags_0(k,j-1,i),5,1)    &
     886                                         )                                     &
     887                          ) * ddy                                              &
     888                        + ( w(k,j,i)   * ( ibit6 + ibit7 + ibit8 )             &
     889                          - w(k-1,j,i) * ( IBITS(wall_flags_0(k-1,j,i),6,1)    &
     890                                         + IBITS(wall_flags_0(k-1,j,i),7,1)    &
     891                                         + IBITS(wall_flags_0(k-1,j,i),8,1)    &
     892                                         )                                     &     
     893                          ) * ddzw(k)
     894
    879895
    880896          tend(k,j,i) = tend(k,j,i) - (                                       &
     
    12471263          DO  k = nzb+1, nzb_max
    12481264
    1249              ibit14 = IBITS(wall_flags_0(k,j,i),14,1)
    1250              ibit13 = IBITS(wall_flags_0(k,j,i),13,1)
    1251              ibit12 = IBITS(wall_flags_0(k,j,i),12,1)
     1265             ibit14 = IBITS(wall_flags_0(k,j-1,i),14,1)
     1266             ibit13 = IBITS(wall_flags_0(k,j-1,i),13,1)
     1267             ibit12 = IBITS(wall_flags_0(k,j-1,i),12,1)
    12521268
    12531269             v_comp      = v(k,j,i) + v(k,j,i-1) - gv
     
    13051321          DO  k = nzb+1, nzb_max
    13061322
    1307              ibit11 = IBITS(wall_flags_0(k,j,i),11,1)
    1308              ibit10 = IBITS(wall_flags_0(k,j,i),10,1)
    1309              ibit9  = IBITS(wall_flags_0(k,j,i),9,1)
     1323             ibit11 = IBITS(wall_flags_0(k,j,i-1),11,1)
     1324             ibit10 = IBITS(wall_flags_0(k,j,i-1),10,1)
     1325             ibit9  = IBITS(wall_flags_0(k,j,i-1),9,1)
    13101326
    13111327             u_comp_l         = u(k,j,i) + u(k,j,i-1) - gu
     
    14811497!--       correction is needed to overcome numerical instabilities introduced
    14821498!--       by a not sufficient reduction of divergences near topography.
    1483           div = ( ( u_comp(k)   - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx       &
    1484                +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy       &
    1485                +  ( w_comp      - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) * ddzw(k)   &
     1499          div = ( ( u_comp(k)       * ( ibit9 + ibit10 + ibit11 )             &
     1500                - ( u(k,j,i)   + u(k,j,i-1)   )                               &
     1501                                    * ( IBITS(wall_flags_0(k,j,i-1),9,1)      &
     1502                                      + IBITS(wall_flags_0(k,j,i-1),10,1)     &
     1503                                      + IBITS(wall_flags_0(k,j,i-1),11,1)     &
     1504                                      )                                       &   
     1505                  ) * ddx                                                     &
     1506               +  ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 )            &
     1507                  - ( v(k,j,i)   + v(k,j,i-1 )  )                             &
     1508                                    * ( IBITS(wall_flags_0(k,j-1,i),12,1)     &
     1509                                      + IBITS(wall_flags_0(k,j-1,i),13,1)     &
     1510                                      + IBITS(wall_flags_0(k,j-1,i),14,1)     &
     1511                                      )                                       &
     1512                  ) * ddy                                                     &
     1513               +  ( w_comp          * ( ibit15 + ibit16 + ibit17 )            &
     1514                - ( w(k-1,j,i) + w(k-1,j,i-1) )                               &
     1515                                    * ( IBITS(wall_flags_0(k-1,j,i),15,1)     &
     1516                                      + IBITS(wall_flags_0(k-1,j,i),16,1)     &
     1517                                      + IBITS(wall_flags_0(k-1,j,i),17,1)     &
     1518                                      )                                       & 
     1519                  ) * ddzw(k)   &
    14861520                ) * 0.5_wp
     1521
    14871522
    14881523          tend(k,j,i) = tend(k,j,i) - (                                       &
     
    17021737          DO  k = nzb+1, nzb_max
    17031738
    1704              ibit20 = IBITS(wall_flags_0(k,j,i),20,1)
    1705              ibit19 = IBITS(wall_flags_0(k,j,i),19,1)
    1706              ibit18 = IBITS(wall_flags_0(k,j,i),18,1)
     1739             ibit20 = IBITS(wall_flags_0(k,j,i-1),20,1)
     1740             ibit19 = IBITS(wall_flags_0(k,j,i-1),19,1)
     1741             ibit18 = IBITS(wall_flags_0(k,j,i-1),18,1)
    17071742
    17081743             u_comp           = u(k,j-1,i) + u(k,j,i) - gu
     
    17601795          DO  k = nzb+1, nzb_max
    17611796
    1762              ibit23 = IBITS(wall_flags_0(k,j,i),23,1)
    1763              ibit22 = IBITS(wall_flags_0(k,j,i),22,1)
    1764              ibit21 = IBITS(wall_flags_0(k,j,i),21,1)
     1797             ibit23 = IBITS(wall_flags_0(k,j-1,i),23,1)
     1798             ibit22 = IBITS(wall_flags_0(k,j-1,i),22,1)
     1799             ibit21 = IBITS(wall_flags_0(k,j-1,i),21,1)
    17651800
    17661801             v_comp_l       = v(k,j,i) + v(k,j-1,i) - gv
     
    19371972!--       correction is needed to overcome numerical instabilities introduced
    19381973!--       by a not sufficient reduction of divergences near topography.
    1939           div = ( ( u_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx       &
    1940                +  ( v_comp(k)   - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy       &
    1941                +  ( w_comp      - ( w(k-1,j-1,i) + w(k-1,j,i) ) ) * ddzw(k)   &
     1974          div = ( ( ( u_comp     + gu )                                       &
     1975                                       * ( ibit18 + ibit19 + ibit20 )         &
     1976                  - ( u(k,j-1,i) + u(k,j,i) )                                 &
     1977                                       * ( IBITS(wall_flags_0(k,j,i-1),18,1)  &
     1978                                         + IBITS(wall_flags_0(k,j,i-1),19,1)  &
     1979                                         + IBITS(wall_flags_0(k,j,i-1),20,1)  &
     1980                                         )                                    &   
     1981                  ) * ddx                                                     &
     1982               +  ( v_comp(k)                                                 &
     1983                                       * ( ibit21 + ibit22 + ibit23 )         &
     1984                - ( v(k,j,i)     + v(k,j-1,i) )                               &
     1985                                       * ( IBITS(wall_flags_0(k,j-1,i),21,1)  &
     1986                                         + IBITS(wall_flags_0(k,j-1,i),22,1)  &
     1987                                         + IBITS(wall_flags_0(k,j-1,i),23,1)  &
     1988                                         )                                    &   
     1989                  ) * ddy                                                     &
     1990               +  ( w_comp                                                    &
     1991                                       * ( ibit24 + ibit25 + ibit26 )         &
     1992                - ( w(k-1,j-1,i) + w(k-1,j,i) )                               &
     1993                                       * ( IBITS(wall_flags_0(k-1,j,i),24,1)  &
     1994                                         + IBITS(wall_flags_0(k-1,j,i),25,1)  &
     1995                                         + IBITS(wall_flags_0(k-1,j,i),26,1)  &
     1996                                         )                                    &
     1997                   ) * ddzw(k)   &
    19421998                ) * 0.5_wp
     1999
    19432000
    19442001          tend(k,j,i) = tend(k,j,i) - (                                       &
     
    21622219
    21632220          DO  k = nzb+1, nzb_max
    2164              ibit32 = IBITS(wall_flags_00(k,j,i),0,1)
    2165              ibit31 = IBITS(wall_flags_0(k,j,i),31,1)
    2166              ibit30 = IBITS(wall_flags_0(k,j,i),30,1)
     2221             ibit32 = IBITS(wall_flags_00(k,j-1,i),0,1)
     2222             ibit31 = IBITS(wall_flags_0(k,j-1,i),31,1)
     2223             ibit30 = IBITS(wall_flags_0(k,j-1,i),30,1)
    21672224
    21682225             v_comp         = v(k+1,j,i) + v(k,j,i) - gv
     
    22202277          DO  k = nzb+1, nzb_max
    22212278
    2222              ibit29 = IBITS(wall_flags_0(k,j,i),29,1)
    2223              ibit28 = IBITS(wall_flags_0(k,j,i),28,1)
    2224              ibit27 = IBITS(wall_flags_0(k,j,i),27,1)
     2279             ibit29 = IBITS(wall_flags_0(k,j,i-1),29,1)
     2280             ibit28 = IBITS(wall_flags_0(k,j,i-1),28,1)
     2281             ibit27 = IBITS(wall_flags_0(k,j,i-1),27,1)
    22252282
    22262283             u_comp           = u(k+1,j,i) + u(k,j,i) - gu
     
    24022459!--       correction is needed to overcome numerical instabilities introduced
    24032460!--       by a not sufficient reduction of divergences near topography.
    2404           div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx         &
    2405               +   ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy         &
    2406               +   ( w_comp      - ( w(k,j,i)   + w(k-1,j,i) ) ) * ddzu(k+1)   &
     2461          div = ( ( ( u_comp + gu ) * ( ibit27 + ibit28 + ibit29 )            &
     2462                  - ( u(k+1,j,i) + u(k,j,i)   )                               &
     2463                                    * ( IBITS(wall_flags_0(k,j,i-1),27,1)     &
     2464                                      + IBITS(wall_flags_0(k,j,i-1),28,1)     &
     2465                                      + IBITS(wall_flags_0(k,j,i-1),29,1)     &
     2466                                      )                                       &
     2467                  ) * ddx                                                     &
     2468              +   ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )            &
     2469                  - ( v(k+1,j,i) + v(k,j,i)   )                               &
     2470                                    * ( IBITS(wall_flags_0(k,j-1,i),30,1)     &
     2471                                      + IBITS(wall_flags_0(k,j-1,i),31,1)     &
     2472                                      + IBITS(wall_flags_00(k,j-1,i),0,1)     &
     2473                                      )                                       &
     2474                  ) * ddy                                                     &
     2475              +   ( w_comp          * ( ibit33 + ibit34 + ibit35 )            &
     2476                - ( w(k,j,i)   + w(k-1,j,i)   )                               &
     2477                                    * ( IBITS(wall_flags_00(k-1,j,i),1,1)     &
     2478                                      + IBITS(wall_flags_00(k-1,j,i),2,1)     &
     2479                                      + IBITS(wall_flags_00(k-1,j,i),3,1)     &
     2480                                      )                                       &
     2481                  ) * ddzu(k+1)   &
    24072482                ) * 0.5_wp
     2483
    24082484
    24092485          tend(k,j,i) = tend(k,j,i) - (                                       &
     
    26312707          DO  k = nzb+1, nzb_max
    26322708
    2633              ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
    2634              ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
    2635              ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
     2709             ibit2 = IBITS(wall_flags_0(k,j,i-1),2,1)
     2710             ibit1 = IBITS(wall_flags_0(k,j,i-1),1,1)
     2711             ibit0 = IBITS(wall_flags_0(k,j,i-1),0,1)
    26362712
    26372713             u_comp                 = u(k,j,i) - u_gtrans
     
    26922768          DO  k = nzb+1, nzb_max
    26932769
    2694              ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
    2695              ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
    2696              ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
     2770             ibit5 = IBITS(wall_flags_0(k,j-1,i),5,1)
     2771             ibit4 = IBITS(wall_flags_0(k,j-1,i),4,1)
     2772             ibit3 = IBITS(wall_flags_0(k,j-1,i),3,1)
    26972773
    26982774             v_comp               = v(k,j,i) - v_gtrans
     
    30063082!--             correction is needed to overcome numerical instabilities caused
    30073083!--             by a not sufficient reduction of divergences near topography.
    3008                 div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx             &
    3009                               + ( v(k,j+1,i) - v(k,j,i)   ) * ddy             &
    3010                               + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
     3084                div   =   ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 )             &
     3085                          - u(k,j,i)   * ( IBITS(wall_flags_0(k,j,i-1),0,1)    &
     3086                                         + IBITS(wall_flags_0(k,j,i-1),1,1)    &
     3087                                         + IBITS(wall_flags_0(k,j,i-1),2,1)    &
     3088                                         )                                     &
     3089                          ) * ddx                                              &
     3090                        + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )             &
     3091                          - v(k,j,i)   * ( IBITS(wall_flags_0(k,j-1,i),3,1)    &
     3092                                         + IBITS(wall_flags_0(k,j-1,i),4,1)    &
     3093                                         + IBITS(wall_flags_0(k,j-1,i),5,1)    &
     3094                                         )                                     &
     3095                          ) * ddy                                              &
     3096                        + ( w(k,j,i)   * ( ibit6 + ibit7 + ibit8 )             &
     3097                          - w(k-1,j,i) * ( IBITS(wall_flags_0(k-1,j,i),6,1)    &
     3098                                         + IBITS(wall_flags_0(k-1,j,i),7,1)    &
     3099                                         + IBITS(wall_flags_0(k-1,j,i),8,1)    &
     3100                                         )                                     &     
     3101                          ) * ddzw(k)
     3102
    30113103
    30123104                tend(k,j,i) = tend(k,j,i) - (                                 &
     
    33883480             DO  k = nzb+1, nzt
    33893481
    3390                 ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
    3391                 ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
    3392                 ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
     3482                ibit2 = IBITS(wall_flags_0(k,j,i-1),2,1)
     3483                ibit1 = IBITS(wall_flags_0(k,j,i-1),1,1)
     3484                ibit0 = IBITS(wall_flags_0(k,j,i-1),0,1)
    33933485
    33943486                u_comp              = u(k,j,i) - u_gtrans
     
    34233515                                                        )
    34243516
     3517                ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
     3518                ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
     3519                ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
     3520
    34253521                u_comp    = u(k,j,i+1) - u_gtrans
    34263522                flux_r    = u_comp * (                                        &
     
    34543550                                             )
    34553551
    3456                 ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
    3457                 ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
    3458                 ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
     3552                ibit5 = IBITS(wall_flags_0(k,j-1,i),5,1)
     3553                ibit4 = IBITS(wall_flags_0(k,j-1,i),4,1)
     3554                ibit3 = IBITS(wall_flags_0(k,j-1,i),3,1)
    34593555
    34603556                v_comp    = v(k,j,i) - v_gtrans
     
    34893585                                             )
    34903586
     3587                ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
     3588                ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
     3589                ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
    34913590
    34923591                v_comp    = v(k,j+1,i) - v_gtrans
     
    37323831!--             correction is needed to overcome numerical instabilities caused
    37333832!--             by a not sufficient reduction of divergences near topography.
    3734                 div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx             &
    3735                               + ( v(k,j+1,i) - v(k,j,i)   ) * ddy             &
    3736                               + ( w(k,j,i)   - w(k-1,j,i) ) * ddzw(k)
     3833                div   =   ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 )             &
     3834                          - u(k,j,i)   * ( IBITS(wall_flags_0(k,j,i-1),0,1)    &
     3835                                         + IBITS(wall_flags_0(k,j,i-1),1,1)    &
     3836                                         + IBITS(wall_flags_0(k,j,i-1),2,1)    &
     3837                                         )                                     &
     3838                          ) * ddx                                              &
     3839                        + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )             &
     3840                          - v(k,j,i)   * ( IBITS(wall_flags_0(k,j-1,i),3,1)    &
     3841                                         + IBITS(wall_flags_0(k,j-1,i),4,1)    &
     3842                                         + IBITS(wall_flags_0(k,j-1,i),5,1)    &
     3843                                         )                                     &
     3844                          ) * ddy                                              &
     3845                        + ( w(k,j,i)   * ( ibit6 + ibit7 + ibit8 )             &
     3846                          - w(k-1,j,i) * ( IBITS(wall_flags_0(k-1,j,i),6,1)    &
     3847                                         + IBITS(wall_flags_0(k-1,j,i),7,1)    &
     3848                                         + IBITS(wall_flags_0(k-1,j,i),8,1)    &
     3849                                         )                                     &     
     3850                          ) * ddzw(k)
     3851
    37373852
    37383853                tend(k,j,i) = - (                                             &
     
    38523967          DO  k = nzb+1, nzb_max
    38533968
    3854              ibit11 = IBITS(wall_flags_0(k,j,i),11,1)
    3855              ibit10 = IBITS(wall_flags_0(k,j,i),10,1)
    3856              ibit9  = IBITS(wall_flags_0(k,j,i),9,1)
     3969             ibit11 = IBITS(wall_flags_0(k,j,i-1),11,1)
     3970             ibit10 = IBITS(wall_flags_0(k,j,i-1),10,1)
     3971             ibit9  = IBITS(wall_flags_0(k,j,i-1),9,1)
    38573972
    38583973             u_comp(k)                = u(k,j,i) + u(k,j,i-1) - gu
     
    39104025          DO  k = nzb+1, nzb_max
    39114026
    3912              ibit14 = IBITS(wall_flags_0(k,j,i),14,1)
    3913              ibit13 = IBITS(wall_flags_0(k,j,i),13,1)
    3914              ibit12 = IBITS(wall_flags_0(k,j,i),12,1)
     4027             ibit14 = IBITS(wall_flags_0(k,j-1,i),14,1)
     4028             ibit13 = IBITS(wall_flags_0(k,j-1,i),13,1)
     4029             ibit12 = IBITS(wall_flags_0(k,j-1,i),12,1)
    39154030
    39164031             v_comp                 = v(k,j,i) + v(k,j,i-1) - gv
     
    40854200!--             correction is needed to overcome numerical instabilities caused
    40864201!--             by a not sufficient reduction of divergences near topography.
    4087                 div = ( ( u_comp(k)   - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx  &
    4088                      +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy  &
    4089                      +  ( w_comp      - ( w(k-1,j,i) + w(k-1,j,i-1) ) )        &
    4090                                                                     * ddzw(k)  &
    4091                       ) * 0.5_wp
     4202                div = ( ( u_comp(k) * ( ibit9 + ibit10 + ibit11 )             &
     4203                - ( u(k,j,i)   + u(k,j,i-1)   )                               &
     4204                                    * ( IBITS(wall_flags_0(k,j,i-1),9,1)      &
     4205                                      + IBITS(wall_flags_0(k,j,i-1),10,1)     &
     4206                                      + IBITS(wall_flags_0(k,j,i-1),11,1)     &
     4207                                      )                                       &   
     4208                  ) * ddx                                                     &
     4209               +  ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 )            &
     4210                  - ( v(k,j,i)   + v(k,j,i-1 )  )                             &
     4211                                    * ( IBITS(wall_flags_0(k,j-1,i),12,1)     &
     4212                                      + IBITS(wall_flags_0(k,j-1,i),13,1)     &
     4213                                      + IBITS(wall_flags_0(k,j-1,i),14,1)     &
     4214                                      )                                       &
     4215                  ) * ddy                                                     &
     4216               +  ( w_comp          * ( ibit15 + ibit16 + ibit17 )            &
     4217                - ( w(k-1,j,i) + w(k-1,j,i-1) )                               &
     4218                                    * ( IBITS(wall_flags_0(k-1,j,i),15,1)     &
     4219                                      + IBITS(wall_flags_0(k-1,j,i),16,1)     &
     4220                                      + IBITS(wall_flags_0(k-1,j,i),17,1)     &
     4221                                      )                                       & 
     4222                  ) * ddzw(k)   &
     4223                ) * 0.5_wp
     4224
     4225
    40924226
    40934227                tend(k,j,i) = tend(k,j,i) - (                                  &
     
    43164450             DO  k = nzb+1, nzt
    43174451
    4318                 ibit11 = IBITS(wall_flags_0(k,j,i),11,1)
    4319                 ibit10 = IBITS(wall_flags_0(k,j,i),10,1)
    4320                 ibit9  = IBITS(wall_flags_0(k,j,i),9,1)
     4452                ibit11 = IBITS(wall_flags_0(k,j,i-1),11,1)
     4453                ibit10 = IBITS(wall_flags_0(k,j,i-1),10,1)
     4454                ibit9  = IBITS(wall_flags_0(k,j,i-1),9,1)
    43214455
    43224456                u_comp_l           = u(k,j,i) + u(k,j,i-1) - gu
     
    43514485                                                         )
    43524486
     4487                ibit11 = IBITS(wall_flags_0(k,j,i),11,1)
     4488                ibit10 = IBITS(wall_flags_0(k,j,i),10,1)
     4489                ibit9  = IBITS(wall_flags_0(k,j,i),9,1)
     4490
    43534491                u_comp    = u(k,j,i+1) + u(k,j,i)
    43544492                flux_r    = ( u_comp   - gu ) * (                           &
     
    43824520                                                     )
    43834521
    4384                 ibit14 = IBITS(wall_flags_0(k,j,i),14,1)
    4385                 ibit13 = IBITS(wall_flags_0(k,j,i),13,1)
    4386                 ibit12 = IBITS(wall_flags_0(k,j,i),12,1)
     4522                ibit14 = IBITS(wall_flags_0(k,j-1,i),14,1)
     4523                ibit13 = IBITS(wall_flags_0(k,j-1,i),13,1)
     4524                ibit12 = IBITS(wall_flags_0(k,j-1,i),12,1)
    43874525
    43884526                v_comp_s                 = v(k,j,i) + v(k,j,i-1) - gv
     
    44184556
    44194557
     4558                ibit14 = IBITS(wall_flags_0(k,j,i),14,1)
     4559                ibit13 = IBITS(wall_flags_0(k,j,i),13,1)
     4560                ibit12 = IBITS(wall_flags_0(k,j,i),12,1)
     4561
    44204562                v_comp    = v(k,j+1,i) + v(k,j+1,i-1) - gv
    44214563                flux_n    = v_comp * (                                       &
     
    45324674!--             correction is needed to overcome numerical instabilities caused
    45334675!--             by a not sufficient reduction of divergences near topography.
    4534                 div = ( ( u_comp      - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx  &
    4535                      +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy  &
    4536                      +  ( w_comp      - ( w(k-1,j,i) + w(k-1,j,i-1) ) )        &
    4537                                                                     * ddzw(k)  &
    4538                       ) * 0.5_wp
     4676                div = ( ( u_comp    * ( ibit9 + ibit10 + ibit11 )             &
     4677                - ( u(k,j,i)   + u(k,j,i-1)   )                               &
     4678                                    * ( IBITS(wall_flags_0(k,j,i-1),9,1)      &
     4679                                      + IBITS(wall_flags_0(k,j,i-1),10,1)     &
     4680                                      + IBITS(wall_flags_0(k,j,i-1),11,1)     &
     4681                                      )                                       &   
     4682                  ) * ddx                                                     &
     4683               +  ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 )            &
     4684                  - ( v(k,j,i)   + v(k,j,i-1 )  )                             &
     4685                                    * ( IBITS(wall_flags_0(k,j-1,i),12,1)     &
     4686                                      + IBITS(wall_flags_0(k,j-1,i),13,1)     &
     4687                                      + IBITS(wall_flags_0(k,j-1,i),14,1)     &
     4688                                      )                                       &
     4689                  ) * ddy                                                     &
     4690               +  ( w_comp          * ( ibit15 + ibit16 + ibit17 )            &
     4691                - ( w(k-1,j,i) + w(k-1,j,i-1) )                               &
     4692                                    * ( IBITS(wall_flags_0(k-1,j,i),15,1)     &
     4693                                      + IBITS(wall_flags_0(k-1,j,i),16,1)     &
     4694                                      + IBITS(wall_flags_0(k-1,j,i),17,1)     &
     4695                                      )                                       & 
     4696                  ) * ddzw(k)   &
     4697                ) * 0.5_wp
     4698
    45394699
    45404700                tend(k,j,i) = - (                                              &
     
    46464806          DO  k = nzb+1, nzb_max
    46474807
    4648              ibit20 = IBITS(wall_flags_0(k,j,i),20,1)
    4649              ibit19 = IBITS(wall_flags_0(k,j,i),19,1)
    4650              ibit18 = IBITS(wall_flags_0(k,j,i),18,1)
     4808             ibit20 = IBITS(wall_flags_0(k,j,i-1),20,1)
     4809             ibit19 = IBITS(wall_flags_0(k,j,i-1),19,1)
     4810             ibit18 = IBITS(wall_flags_0(k,j,i-1),18,1)
    46514811
    46524812             u_comp                   = u(k,j-1,i) + u(k,j,i) - gu
     
    47044864          DO  k = nzb+1, nzb_max
    47054865
    4706              ibit23 = IBITS(wall_flags_0(k,j,i),23,1)
    4707              ibit22 = IBITS(wall_flags_0(k,j,i),22,1)
    4708              ibit21 = IBITS(wall_flags_0(k,j,i),21,1)
     4866             ibit23 = IBITS(wall_flags_0(k,j-1,i),23,1)
     4867             ibit22 = IBITS(wall_flags_0(k,j-1,i),22,1)
     4868             ibit21 = IBITS(wall_flags_0(k,j-1,i),21,1)
    47094869
    47104870             v_comp(k)              = v(k,j,i) + v(k,j-1,i) - gv
     
    48785038!--             correction is needed to overcome numerical instabilities caused
    48795039!--             by a not sufficient reduction of divergences near topography.
    4880                 div = ( ( u_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx &
    4881                      +  ( v_comp(k)   - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy &
    4882                      +  ( w_comp      - ( w(k-1,j-1,i) + w(k-1,j,i) )         &
    4883                                                                   ) * ddzw(k) &
    4884                       ) * 0.5_wp
     5040                div = ( ( ( u_comp     + gu )                                 &
     5041                                       * ( ibit18 + ibit19 + ibit20 )         &
     5042                - ( u(k,j-1,i)   + u(k,j,i) )                                 &
     5043                                       * ( IBITS(wall_flags_0(k,j,i-1),18,1)  &
     5044                                         + IBITS(wall_flags_0(k,j,i-1),19,1)  &
     5045                                         + IBITS(wall_flags_0(k,j,i-1),20,1)  &
     5046                                         )                                    &   
     5047                  ) * ddx                                                     &
     5048               +  ( v_comp(k)                                                 &
     5049                                       * ( ibit21 + ibit22 + ibit23 )         &
     5050                - ( v(k,j,i)     + v(k,j-1,i) )                               &
     5051                                       * ( IBITS(wall_flags_0(k,j-1,i),21,1)  &
     5052                                         + IBITS(wall_flags_0(k,j-1,i),22,1)  &
     5053                                         + IBITS(wall_flags_0(k,j-1,i),23,1)  &
     5054                                         )                                    &   
     5055                  ) * ddy                                                     &
     5056               +  ( w_comp                                                    &
     5057                                       * ( ibit24 + ibit25 + ibit26 )         &
     5058                - ( w(k-1,j-1,i) + w(k-1,j,i) )                               &
     5059                                       * ( IBITS(wall_flags_0(k-1,j,i),24,1)  &
     5060                                         + IBITS(wall_flags_0(k-1,j,i),25,1)  &
     5061                                         + IBITS(wall_flags_0(k-1,j,i),26,1)  &
     5062                                         )                                    &
     5063                   ) * ddzw(k)   &
     5064                ) * 0.5_wp
     5065
    48855066
    48865067                tend(k,j,i) = tend(k,j,i) - (                                 &
     
    51205301             DO  k = nzb+1, nzt
    51215302
    5122                 ibit20 = IBITS(wall_flags_0(k,j,i),20,1)
    5123                 ibit19 = IBITS(wall_flags_0(k,j,i),19,1)
    5124                 ibit18 = IBITS(wall_flags_0(k,j,i),18,1)
     5303                ibit20 = IBITS(wall_flags_0(k,j,i-1),20,1)
     5304                ibit19 = IBITS(wall_flags_0(k,j,i-1),19,1)
     5305                ibit18 = IBITS(wall_flags_0(k,j,i-1),18,1)
    51255306
    51265307                u_comp_l                 = u(k,j-1,i) + u(k,j,i) - gu
     
    51555336                                                           )
    51565337
     5338                ibit20 = IBITS(wall_flags_0(k,j,i),20,1)
     5339                ibit19 = IBITS(wall_flags_0(k,j,i),19,1)
     5340                ibit18 = IBITS(wall_flags_0(k,j,i),18,1)
     5341
    51575342                u_comp    = u(k,j-1,i+1) + u(k,j,i+1) - gu
    51585343                flux_r    = u_comp * (                                       &
     
    51865371                                              )
    51875372
    5188                 ibit23 = IBITS(wall_flags_0(k,j,i),23,1)
    5189                 ibit22 = IBITS(wall_flags_0(k,j,i),22,1)
    5190                 ibit21 = IBITS(wall_flags_0(k,j,i),21,1)
     5373                ibit23 = IBITS(wall_flags_0(k,j-1,i),23,1)
     5374                ibit22 = IBITS(wall_flags_0(k,j-1,i),22,1)
     5375                ibit21 = IBITS(wall_flags_0(k,j-1,i),21,1)
    51915376
    51925377
     
    52225407                                                          )
    52235408
     5409                ibit23 = IBITS(wall_flags_0(k,j,i),23,1)
     5410                ibit22 = IBITS(wall_flags_0(k,j,i),22,1)
     5411                ibit21 = IBITS(wall_flags_0(k,j,i),21,1)
     5412
    52245413                v_comp = v(k,j+1,i) + v(k,j,i)
    52255414                flux_n = ( v_comp - gv ) * (                                 &
     
    53365525!--             correction is needed to overcome numerical instabilities caused
    53375526!--             by a not sufficient reduction of divergences near topography.
    5338                 div = ( ( u_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx &
    5339                      +  ( v_comp      - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy &
    5340                      +  ( w_comp      - ( w(k-1,j-1,i) + w(k-1,j,i) )         &
    5341                                                                   ) * ddzw(k) &
    5342                       ) * 0.5_wp
     5527                div = ( ( ( u_comp     + gu )                                 &
     5528                                       * ( ibit18 + ibit19 + ibit20 )         &
     5529                - ( u(k,j-1,i)   + u(k,j,i) )                                 &
     5530                                       * ( IBITS(wall_flags_0(k,j,i-1),18,1)  &
     5531                                         + IBITS(wall_flags_0(k,j,i-1),19,1)  &
     5532                                         + IBITS(wall_flags_0(k,j,i-1),20,1)  &
     5533                                         )                                    &   
     5534                  ) * ddx                                                     &
     5535               +  ( v_comp                                                    &
     5536                                       * ( ibit21 + ibit22 + ibit23 )         &
     5537                - ( v(k,j,i)     + v(k,j-1,i) )                               &
     5538                                       * ( IBITS(wall_flags_0(k,j-1,i),21,1)  &
     5539                                         + IBITS(wall_flags_0(k,j-1,i),22,1)  &
     5540                                         + IBITS(wall_flags_0(k,j-1,i),23,1)  &
     5541                                         )                                    &   
     5542                  ) * ddy                                                     &
     5543               +  ( w_comp                                                    &
     5544                                       * ( ibit24 + ibit25 + ibit26 )         &
     5545                - ( w(k-1,j-1,i) + w(k-1,j,i) )                               &
     5546                                       * ( IBITS(wall_flags_0(k-1,j,i),24,1)  &
     5547                                         + IBITS(wall_flags_0(k-1,j,i),25,1)  &
     5548                                         + IBITS(wall_flags_0(k-1,j,i),26,1)  &
     5549                                         )                                    &
     5550                   ) * ddzw(k)   &
     5551                ) * 0.5_wp
     5552
    53435553
    53445554                tend(k,j,i) = - (                                              &
     
    54525662          DO  k = nzb+1, nzb_max
    54535663
    5454              ibit29 = IBITS(wall_flags_0(k,j,i),29,1)
    5455              ibit28 = IBITS(wall_flags_0(k,j,i),28,1)
    5456              ibit27 = IBITS(wall_flags_0(k,j,i),27,1)
     5664             ibit29 = IBITS(wall_flags_0(k,j,i-1),29,1)
     5665             ibit28 = IBITS(wall_flags_0(k,j,i-1),28,1)
     5666             ibit27 = IBITS(wall_flags_0(k,j,i-1),27,1)
    54575667
    54585668             u_comp                   = u(k+1,j,i) + u(k,j,i) - gu
     
    55105720          DO  k = nzb+1, nzb_max
    55115721
    5512              ibit32 = IBITS(wall_flags_00(k,j,i),0,1)
    5513              ibit31 = IBITS(wall_flags_0(k,j,i),31,1)
    5514              ibit30 = IBITS(wall_flags_0(k,j,i),30,1)
     5722             ibit32 = IBITS(wall_flags_00(k,j-1,i),0,1)
     5723             ibit31 = IBITS(wall_flags_0(k,j-1,i),31,1)
     5724             ibit30 = IBITS(wall_flags_0(k,j-1,i),30,1)
    55155725
    55165726             v_comp                 = v(k+1,j,i) + v(k,j,i) - gv
     
    56905900!--             correction is needed to overcome numerical instabilities caused
    56915901!--             by a not sufficient reduction of divergences near topography.
    5692                 div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx  &
    5693                     +   ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy  &
    5694                     +   ( w_comp      - ( w(k,j,i)   + w(k-1,j,i) ) )        &
    5695                                                                  * ddzu(k+1) &
    5696                       ) * 0.5_wp
     5902                div = ( ( ( u_comp + gu ) * ( ibit27 + ibit28 + ibit29 )      &
     5903                  - ( u(k+1,j,i) + u(k,j,i)   )                               &
     5904                                    * ( IBITS(wall_flags_0(k,j,i-1),27,1)     &
     5905                                      + IBITS(wall_flags_0(k,j,i-1),28,1)     &
     5906                                      + IBITS(wall_flags_0(k,j,i-1),29,1)     &
     5907                                      )                                       &
     5908                  ) * ddx                                                     &
     5909              +   ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )            &
     5910                  - ( v(k+1,j,i) + v(k,j,i)   )                               &
     5911                                    * ( IBITS(wall_flags_0(k,j-1,i),30,1)     &
     5912                                      + IBITS(wall_flags_0(k,j-1,i),31,1)     &
     5913                                      + IBITS(wall_flags_00(k,j-1,i),0,1)     &
     5914                                      )                                       &
     5915                  ) * ddy                                                     &
     5916              +   ( w_comp          * ( ibit33 + ibit34 + ibit35 )            &
     5917                - ( w(k,j,i)   + w(k-1,j,i)   )                               &
     5918                                    * ( IBITS(wall_flags_00(k-1,j,i),1,1)     &
     5919                                      + IBITS(wall_flags_00(k-1,j,i),2,1)     &
     5920                                      + IBITS(wall_flags_00(k-1,j,i),3,1)     &
     5921                                      )                                       &
     5922                  ) * ddzu(k+1)   &
     5923                ) * 0.5_wp
     5924
     5925
    56975926
    56985927                tend(k,j,i) = tend(k,j,i) - (                                 &
     
    59036132             DO  k = nzb+1, nzt
    59046133
    5905                 ibit27 = IBITS(wall_flags_0(k,j,i),27,1)
    5906                 ibit28 = IBITS(wall_flags_0(k,j,i),28,1)
    5907                 ibit29 = IBITS(wall_flags_0(k,j,i),29,1)
     6134                ibit27 = IBITS(wall_flags_0(k,j,i-1),27,1)
     6135                ibit28 = IBITS(wall_flags_0(k,j,i-1),28,1)
     6136                ibit29 = IBITS(wall_flags_0(k,j,i-1),29,1)
    59086137
    59096138                u_comp_l                 = u(k+1,j,i) + u(k,j,i) - gu
     
    59386167                                                            )
    59396168
     6169                ibit27 = IBITS(wall_flags_0(k,j,i),27,1)
     6170                ibit28 = IBITS(wall_flags_0(k,j,i),28,1)
     6171                ibit29 = IBITS(wall_flags_0(k,j,i),29,1)
     6172
    59406173                u_comp    = u(k+1,j,i+1) + u(k,j,i+1) - gu
    59416174                flux_r    = u_comp * (                                       &
     
    59686201                                 ( w(k,j,i+3) - w(k,j,i-2) )                 &
    59696202                                              )
    5970                 ibit32 = IBITS(wall_flags_00(k,j,i),0,1)
    5971                 ibit31 = IBITS(wall_flags_0(k,j,i),31,1)
    5972                 ibit30 = IBITS(wall_flags_0(k,j,i),30,1)
     6203                ibit32 = IBITS(wall_flags_00(k,j-1,i),0,1)
     6204                ibit31 = IBITS(wall_flags_0(k,j-1,i),31,1)
     6205                ibit30 = IBITS(wall_flags_0(k,j-1,i),30,1)
    59736206
    59746207                v_comp_s               = v(k+1,j,i) + v(k,j,i) - gv
     
    60036236                                                        )
    60046237
     6238                ibit32 = IBITS(wall_flags_00(k,j,i),0,1)
     6239                ibit31 = IBITS(wall_flags_0(k,j,i),31,1)
     6240                ibit30 = IBITS(wall_flags_0(k,j,i),30,1)
     6241
    60056242                v_comp    = v(k+1,j+1,i) + v(k,j+1,i) - gv
    60066243                flux_n    = v_comp * (                                       &
     
    61186355!--             correction is needed to overcome numerical instabilities caused
    61196356!--             by a not sufficient reduction of divergences near topography.
    6120                 div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx  &
    6121                     +   ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy  &
    6122                     +   ( w_comp      - ( w(k,j,i)   + w(k-1,j,i) ) )        &
    6123                                                                  * ddzu(k+1) &
    6124                       ) * 0.5_wp
     6357                div = ( ( ( u_comp + gu ) * ( ibit27 + ibit28 + ibit29 )      &
     6358                  - ( u(k+1,j,i) + u(k,j,i)   )                               &
     6359                                    * ( IBITS(wall_flags_0(k,j,i-1),27,1)     &
     6360                                      + IBITS(wall_flags_0(k,j,i-1),28,1)     &
     6361                                      + IBITS(wall_flags_0(k,j,i-1),29,1)     &
     6362                                      )                                       &
     6363                  ) * ddx                                                     &
     6364              +   ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )            &
     6365                  - ( v(k+1,j,i) + v(k,j,i)   )                               &
     6366                                    * ( IBITS(wall_flags_0(k,j-1,i),30,1)     &
     6367                                      + IBITS(wall_flags_0(k,j-1,i),31,1)     &
     6368                                      + IBITS(wall_flags_00(k,j-1,i),0,1)     &
     6369                                      )                                       &
     6370                  ) * ddy                                                     &
     6371              +   ( w_comp          * ( ibit33 + ibit34 + ibit35 )            &
     6372                - ( w(k,j,i)   + w(k-1,j,i)   )                               &
     6373                                    * ( IBITS(wall_flags_00(k-1,j,i),1,1)     &
     6374                                      + IBITS(wall_flags_00(k-1,j,i),2,1)     &
     6375                                      + IBITS(wall_flags_00(k-1,j,i),3,1)     &
     6376                                      )                                       &
     6377                  ) * ddzu(k+1)   &
     6378                ) * 0.5_wp
     6379
    61256380
    61266381                tend(k,j,i) = - (                                                &
Note: See TracChangeset for help on using the changeset viewer.