Ignore:
Timestamp:
Aug 3, 2017 2:24:56 PM (4 years ago)
Author:
knoop
Message:

Bugfix for topography usage with anelastic approximation and boussinesq approximation with air density != 1

File:
1 edited

Legend:

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

    r2292 r2329  
    2525! -----------------
    2626! $Id$
     27! Bugfix concerning density in divergence correction close to buildings
     28!
     29! 2292 2017-06-20 09:51:42Z schwenkel
    2730! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
    2831! includes two more prognostic equations for cloud drop concentration (nc) 
     
    10941097
    10951098       USE arrays_3d,                                                          &
    1096            ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
     1099           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air_zw
    10971100
    10981101       USE constants,                                                          &
     
    14221425                                         + IBITS(advc_flags_1(k,j,i-1),2,1)    &
    14231426                                         )                                     &
    1424                           ) * rho_air(k) * ddx                                 &
     1427                          ) * ddx                                              &
    14251428                        + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )             &
    14261429                          - v(k,j,i)   * ( IBITS(advc_flags_1(k,j-1,i),3,1)    &
     
    14281431                                         + IBITS(advc_flags_1(k,j-1,i),5,1)    &
    14291432                                         )                                     &
    1430                           ) * rho_air(k) * ddy                                 &
     1433                          ) * ddy                                              &
    14311434                        + ( w(k,j,i) * rho_air_zw(k) *                         &
    14321435                                         ( ibit6 + ibit7 + ibit8 )             &
     
    14361439                                         + IBITS(advc_flags_1(k-1,j,i),8,1)    &
    14371440                                         )                                     &     
    1438                           ) * ddzw(k)
     1441                          ) * drho_air(k) * ddzw(k)
    14391442
    14401443
     
    15251528!--       correction is needed to overcome numerical instabilities introduced
    15261529!--       by a not sufficient reduction of divergences near topography.
    1527           div         =   ( u(k,j,i+1) - u(k,j,i)   ) * rho_air(k) * ddx      &
    1528                         + ( v(k,j+1,i) - v(k,j,i)   ) * rho_air(k) * ddy      &
     1530          div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx                   &
     1531                        + ( v(k,j+1,i) - v(k,j,i)   ) * ddy                   &
    15291532                        + ( w(k,j,i)   * rho_air_zw(k) -                      &
    1530                             w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k)
     1533                            w(k-1,j,i) * rho_air_zw(k-1)                      &
     1534                          ) * drho_air(k) * ddzw(k)
    15311535
    15321536          tend(k,j,i) = tend(k,j,i) - (                                       &
     
    16591663       USE arrays_3d,                                                         &
    16601664           ONLY:  ddzw, diss_l_u, diss_s_u, flux_l_u, flux_s_u, tend, u, v, w,&
    1661                   drho_air, rho_air, rho_air_zw
     1665                  drho_air, rho_air_zw
    16621666
    16631667       USE constants,                                                         &
     
    19631967                                      + IBITS(advc_flags_1(k,j,i-1),11,1)     &
    19641968                                      )                                       &
    1965                   ) * rho_air(k) * ddx                                        &
     1969                  ) * ddx                                                     &
    19661970               +  ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 )            &
    19671971                  - ( v(k,j,i)   + v(k,j,i-1 )  )                             &
     
    19701974                                      + IBITS(advc_flags_1(k,j-1,i),14,1)     &
    19711975                                      )                                       &
    1972                   ) * rho_air(k) * ddy                                        &
     1976                  ) * ddy                                                     &
    19731977               +  ( w_comp * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 )     &
    19741978                - ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1)             &
     
    19771981                                      + IBITS(advc_flags_1(k-1,j,i),17,1)     &
    19781982                                      )                                       & 
    1979                   ) * ddzw(k)   &
     1983                  ) * drho_air(k) * ddzw(k)                                   &
    19801984                ) * 0.5_wp
    19811985
     
    20872091!--       by a not sufficient reduction of divergences near topography.
    20882092          div = ( ( u_comp(k)   - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx       &
    2089                                                                   * rho_air(k)&
    20902093               +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy       &
    2091                                                                   * rho_air(k)&
    20922094               +  (   w_comp                      * rho_air_zw(k)   -         &
    20932095                    ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1)           &
    2094                   ) * ddzw(k)                                                 &
     2096                  ) * drho_air(k) * ddzw(k)                                   &
    20952097                ) * 0.5_wp
    20962098
     
    21512153       USE arrays_3d,                                                          &
    21522154           ONLY:  ddzw, diss_l_v, diss_s_v, flux_l_v, flux_s_v, tend, u, v, w, &
    2153                   drho_air, rho_air, rho_air_zw
     2155                  drho_air, rho_air_zw
    21542156
    21552157       USE constants,                                                          &
     
    24582460                                         + IBITS(advc_flags_1(k,j,i-1),20,1)  &
    24592461                                         )                                    &
    2460                   ) * rho_air(k) * ddx                                        &
     2462                  ) * ddx                                                     &
    24612463               +  ( v_comp(k)                                                 &
    24622464                                       * ( ibit21 + ibit22 + ibit23 )         &
     
    24662468                                         + IBITS(advc_flags_1(k,j-1,i),23,1)  &
    24672469                                         )                                    &
    2468                   ) * rho_air(k) * ddy                                        &
     2470                  ) * ddy                                                     &
    24692471               +  ( w_comp * rho_air_zw(k) * ( ibit24 + ibit25 + ibit26 )     &
    24702472                - ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)             &
     
    24732475                                         + IBITS(advc_flags_1(k-1,j,i),26,1)  &
    24742476                                         )                                    &
    2475                    ) * ddzw(k)   &
     2477                   ) * drho_air(k) * ddzw(k)                                  &
    24762478                ) * 0.5_wp
    24772479
     
    25882590!--       by a not sufficient reduction of divergences near topography.
    25892591          div = ( ( u_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx       &
    2590                                                                   * rho_air(k)&
    25912592               +  ( v_comp(k)   - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy       &
    2592                                                                   * rho_air(k)&
    25932593               +  (   w_comp                      * rho_air_zw(k)   -         &
    25942594                    ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)           &
    2595                   ) * ddzw(k)                                                 &
     2595                  ) * drho_air(k) * ddzw(k)                                   &
    25962596                ) * 0.5_wp
    25972597
     
    26522652       USE arrays_3d,                                                         &
    26532653           ONLY:  ddzu, diss_l_w, diss_s_w, flux_l_w, flux_s_w, tend, u, v, w,&
    2654                   drho_air_zw, rho_air, rho_air_zw
     2654                  drho_air_zw, rho_air
    26552655
    26562656       USE constants,                                                         &
     
    29622962                                      + IBITS(advc_flags_1(k,j,i-1),29,1)     &
    29632963                                      )                                       &
    2964                   ) * rho_air_zw(k) * ddx                                     &
     2964                  ) * ddx                                                     &
    29652965              +   ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )            &
    29662966                  - ( v(k+1,j,i) + v(k,j,i)   )                               &
     
    29692969                                      + IBITS(advc_flags_2(k,j-1,i),0,1)      &
    29702970                                      )                                       &
    2971                   ) * rho_air_zw(k) * ddy                                     &
     2971                  ) * ddy                                                     &
    29722972              +   ( w_comp * rho_air(k+1) * ( ibit33 + ibit34 + ibit35 )      &
    29732973                - ( w(k,j,i)   + w(k-1,j,i)   ) * rho_air(k)                  &
     
    29762976                                      + IBITS(advc_flags_2(k-1,j,i),3,1)      &
    29772977                                      )                                       &
    2978                   ) * ddzu(k+1)   &
     2978                  ) * drho_air_zw(k) * ddzu(k+1)                              &
    29792979                ) * 0.5_wp
    29802980
     
    30783078!--       by a not sufficient reduction of divergences near topography.
    30793079          div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx         &
    3080                                                               * rho_air_zw(k) &
    30813080              +   ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy         &
    3082                                                               * rho_air_zw(k) &
    30833081              +   (   w_comp                    * rho_air(k+1) -              &
    30843082                    ( w(k,j,i)   + w(k-1,j,i) ) * rho_air(k)                  &
    3085                   ) * ddzu(k+1)                                               &
     3083                  ) * drho_air_zw(k) * ddzu(k+1)                              &
    30863084                ) * 0.5_wp
    30873085
     
    31273125
    31283126       USE arrays_3d,                                                         &
    3129            ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
     3127           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air_zw
    31303128
    31313129       USE constants,                                                         &
     
    34483446                                         + IBITS(advc_flags_1(k,j,i-1),2,1)    &
    34493447                                         )                                     &
    3450                           ) * rho_air(k) * ddx                                 &
     3448                          ) * ddx                                              &
    34513449                        + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 )             &
    34523450                          - v(k,j,i)   * ( IBITS(advc_flags_1(k,j-1,i),3,1)    &
     
    34543452                                         + IBITS(advc_flags_1(k,j-1,i),5,1)    &
    34553453                                         )                                     &
    3456                           ) * rho_air(k) * ddy                                 &
     3454                          ) * ddy                                              &
    34573455                        + ( w(k,j,i) * rho_air_zw(k) *                         &
    34583456                                         ( ibit6 + ibit7 + ibit8 )             &
     
    34623460                                         + IBITS(advc_flags_1(k-1,j,i),8,1)    &
    34633461                                         )                                     &     
    3464                           ) * ddzw(k)
     3462                          ) * drho_air(k) * ddzw(k)
    34653463
    34663464
     
    35493547!--             correction is needed to overcome numerical instabilities introduced
    35503548!--             by a not sufficient reduction of divergences near topography.
    3551                 div         =   ( u(k,j,i+1) - u(k,j,i)   ) * rho_air(k) * ddx &
    3552                               + ( v(k,j+1,i) - v(k,j,i)   ) * rho_air(k) * ddy &
     3549                div         =   ( u(k,j,i+1) - u(k,j,i)   ) * ddx              &
     3550                              + ( v(k,j+1,i) - v(k,j,i)   ) * ddy              &
    35533551                              + ( w(k,j,i)   * rho_air_zw(k) -                 &
    3554                                   w(k-1,j,i) * rho_air_zw(k-1) ) * ddzw(k)
     3552                                  w(k-1,j,i) * rho_air_zw(k-1)                 &
     3553                                ) * drho_air(k) * ddzw(k)
    35553554
    35563555                tend(k,j,i) = tend(k,j,i) - (                                 &
     
    36823681
    36833682       USE arrays_3d,                                                          &
    3684            ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
     3683           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air_zw
    36853684
    36863685       USE constants,                                                          &
     
    39903989                                      + IBITS(advc_flags_1(k,j,i-1),11,1)     &
    39913990                                      )                                       &
    3992                   ) * rho_air(k) * ddx                                        &
     3991                  ) * ddx                                                     &
    39933992               +  ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 )            &
    39943993                  - ( v(k,j,i)   + v(k,j,i-1 )  )                             &
     
    39973996                                      + IBITS(advc_flags_1(k,j-1,i),14,1)     &
    39983997                                      )                                       &
    3999                   ) * rho_air(k) * ddy                                        &
     3998                  ) * ddy                                                     &
    40003999               +  ( w_comp * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 )     &
    40014000                - ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1)             &
     
    40044003                                      + IBITS(advc_flags_1(k-1,j,i),17,1)     &
    40054004                                      )                                       & 
    4006                   ) * ddzw(k)   &
     4005                  ) * drho_air(k) * ddzw(k)                                   &
    40074006                ) * 0.5_wp
    40084007
     
    41164115!--             by a not sufficient reduction of divergences near topography.
    41174116                div = ( ( u_comp(k)   - ( u(k,j,i)   + u(k,j,i-1)   ) ) * ddx &
    4118                                                                   * rho_air(k)&
    41194117                     +  ( v_comp + gv - ( v(k,j,i)   + v(k,j,i-1 )  ) ) * ddy &
    4120                                                                   * rho_air(k)&
    41214118                     +  (   w_comp                      * rho_air_zw(k) -     &
    41224119                          ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1)     &
    4123                         ) * ddzw(k)                                           &
     4120                        ) * drho_air(k) * ddzw(k)                             &
    41244121                      ) * 0.5_wp
    41254122
     
    41784175
    41794176       USE arrays_3d,                                                          &
    4180            ONLY:  ddzw, drho_air, tend, u, v, w, rho_air, rho_air_zw
     4177           ONLY:  ddzw, drho_air, tend, u, v, w, rho_air_zw
    41814178
    41824179       USE constants,                                                          &
     
    44864483                                         + IBITS(advc_flags_1(k,j,i-1),20,1)  &
    44874484                                         )                                    &
    4488                   ) * rho_air(k) * ddx                                        &
     4485                  ) * ddx                                                     &
    44894486               +  ( v_comp(k)                                                 &
    44904487                                       * ( ibit21 + ibit22 + ibit23 )         &
     
    44944491                                         + IBITS(advc_flags_1(k,j-1,i),23,1)  &
    44954492                                         )                                    &
    4496                   ) * rho_air(k) * ddy                                        &
     4493                  ) * ddy                                                     &
    44974494               +  ( w_comp * rho_air_zw(k)                                    &
    44984495                                       * ( ibit24 + ibit25 + ibit26 )         &
     
    45024499                                         + IBITS(advc_flags_1(k-1,j,i),26,1)  &
    45034500                                         )                                    &
    4504                    ) * ddzw(k)   &
     4501                   ) * drho_air(k) * ddzw(k)                                  &
    45054502                ) * 0.5_wp
    45064503
     
    46194616!--             by a not sufficient reduction of divergences near topography.
    46204617                div = ( ( u_comp + gu - ( u(k,j-1,i)   + u(k,j,i)   ) ) * ddx &
    4621                                                                   * rho_air(k)&
    46224618                     +  ( v_comp(k)   - ( v(k,j,i)     + v(k,j-1,i) ) ) * ddy &
    4623                                                                   * rho_air(k)&
    46244619                     +  (   w_comp                      * rho_air_zw(k) -     &
    46254620                          ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1)     &
    4626                         ) * ddzw(k)                                           &
     4621                        ) * drho_air(k) * ddzw(k)                             &
    46274622                      ) * 0.5_wp
    46284623 
     
    46854680
    46864681       USE arrays_3d,                                                          &
    4687            ONLY:  ddzu, drho_air_zw, tend, u, v, w, rho_air, rho_air_zw
     4682           ONLY:  ddzu, drho_air_zw, tend, u, v, w, rho_air
    46884683
    46894684       USE constants,                                                          &
     
    49984993                                      + IBITS(advc_flags_1(k,j,i-1),29,1)     &
    49994994                                      )                                       &
    5000                   ) * rho_air_zw(k) * ddx                                     &
     4995                  ) * ddx                                                     &
    50014996              +   ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 )            &
    50024997                  - ( v(k+1,j,i) + v(k,j,i)   )                               &
     
    50055000                                      + IBITS(advc_flags_2(k,j-1,i),0,1)      &
    50065001                                      )                                       &
    5007                   ) * rho_air_zw(k) * ddy                                     &
     5002                  ) * ddy                                                     &
    50085003              +   ( w_comp * rho_air(k+1) * ( ibit33 + ibit34 + ibit35 )      &
    50095004                - ( w(k,j,i)   + w(k-1,j,i)   ) * rho_air(k)                  &
     
    50125007                                      + IBITS(advc_flags_2(k-1,j,i),3,1)      &
    50135008                                      )                                       &
    5014                   ) * ddzu(k+1)   &
     5009                  ) * drho_air_zw(k) * ddzu(k+1)                              &
    50155010                ) * 0.5_wp
    50165011
     
    51165111!--             by a not sufficient reduction of divergences near topography.
    51175112                div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i)   ) ) * ddx  &
    5118                                                              * rho_air_zw(k) &
    51195113                    +   ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i)   ) ) * ddy  &
    5120                                                              * rho_air_zw(k) &
    51215114                    +   (   w_comp                    * rho_air(k+1) -       &
    51225115                          ( w(k,j,i)   + w(k-1,j,i) ) * rho_air(k)           &
    5123                         ) * ddzu(k+1)                                        &
     5116                        ) * drho_air_zw(k) * ddzu(k+1)                       &
    51245117                      ) * 0.5_wp
    51255118
Note: See TracChangeset for help on using the changeset viewer.