Changeset 2329 for palm/trunk/SOURCE


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

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

Location:
palm/trunk/SOURCE
Files:
4 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
  • palm/trunk/SOURCE/check_parameters.f90

    r2320 r2329  
    2525! -----------------
    2626! $Id$
     27! Bugfix: index corrected for rho_air and rho_air_zw output
     28!
     29! 2320 2017-07-21 12:47:43Z suehring
    2730! Modularize large-scale forcing and nudging
    2831!
     
    27482751
    27492752          CASE ( 'rho_air' )
    2750              dopr_index(i)  = 121
     2753             dopr_index(i)  = 119
    27512754             dopr_unit(i)   = 'kg/m3'
    2752              hom(:,2,121,:) = SPREAD( zu, 2, statistic_regions+1 )
     2755             hom(:,2,119,:) = SPREAD( zu, 2, statistic_regions+1 )
    27532756
    27542757          CASE ( 'rho_air_zw' )
    2755              dopr_index(i)  = 122
     2758             dopr_index(i)  = 120
    27562759             dopr_unit(i)   = 'kg/m3'
    2757              hom(:,2,122,:) = SPREAD( zw, 2, statistic_regions+1 )
     2760             hom(:,2,120,:) = SPREAD( zw, 2, statistic_regions+1 )
    27582761
    27592762          CASE ( 'nc' )
  • palm/trunk/SOURCE/init_3d_model.f90

    r2327 r2329  
    2525! -----------------
    2626! $Id$
     27! Removed temporary bugfix (r2327) as bug is properly resolved by this revision
     28!
     29! 2327 2017-08-02 07:40:57Z maronga
    2730! Temporary bugfix
    2831!
     
    346349!> or
    347350!> c) read values of a previous run
    348 !> @todo fix crashes for runs with topography and density /= 1
    349351!------------------------------------------------------------------------------!
    350352 SUBROUTINE init_3d_model
     
    737739       rho_air_zw(nzt+1)  = rho_air_zw(nzt)                                    &
    738740                            + 2.0_wp * ( rho_air(nzt+1) - rho_air_zw(nzt)  )
    739 
    740 !
    741 !--    Temporary workaround as runs with topography do not work with density /= 1
    742        rho_air    = 1.0_wp
    743        rho_air_zw = 1.0_wp
    744 
    745     ENDIF
    746 
    747 
    748 
     741    ENDIF
    749742
    750743!-- compute the inverse density array in order to avoid expencive divisions
  • palm/trunk/SOURCE/production_e.f90

    r2233 r2329  
    2525! -----------------
    2626! $Id$
     27! Bugfix: added division by density as kinematic fluxes are needed
     28!
     29! 2233 2017-05-30 18:08:54Z suehring
    2730!
    2831! 2232 2017-05-30 17:47:52Z suehring
     
    135138
    136139       USE arrays_3d,                                                          &
    137            ONLY:  ddzw, dd2zu, kh, km, prho, pt, q, ql, tend, u, v, vpt, w
     140           ONLY:  ddzw, dd2zu, drho_air_zw, kh, km, prho, pt, q, ql, tend, u,  &
     141                  v, vpt, w
    138142
    139143       USE cloud_parameters,                                                   &
     
    531535                                  k = surf_def_h(l)%k(m)
    532536                                  tend(k,j,i) = tend(k,j,i) + g / pt_reference &
     537                                                         * drho_air_zw(k-1)    &
    533538                                                         * surf_def_h(l)%shf(m)   
    534539                               ENDDO   
     
    541546                               k = surf_lsm_h%k(m)
    542547                               tend(k,j,i) = tend(k,j,i) + g / pt_reference    &
     548                                                         * drho_air_zw(k-1)    &
    543549                                                         * surf_lsm_h%shf(m)   
    544550                            ENDDO
     
    550556                               k = surf_usm_h%k(m)
    551557                               tend(k,j,i) = tend(k,j,i) + g / pt_reference    &
     558                                                         * drho_air_zw(k-1)    &
    552559                                                         * surf_usm_h%shf(m)   
    553560                            ENDDO                         
     
    559566                            DO  m = surf_s, surf_e
    560567                               k = surf_def_h(2)%k(m)
    561                                tend(k,j,i) = tend(k,j,i) + g / pt_reference *     &
    562                                                            surf_def_h(2)%shf(m)
     568                               tend(k,j,i) = tend(k,j,i) + g / pt_reference    &
     569                                                         * drho_air_zw(k-1)    &
     570                                                         * surf_def_h(2)%shf(m)
    563571                            ENDDO
    564572                         ENDIF
     
    613621                                  k = surf_def_h(l)%k(m)
    614622                                  tend(k,j,i) = tend(k,j,i) + g / pt_reference &
     623                                                         * drho_air_zw(k-1)    &
    615624                                                         * surf_def_h(l)%shf(m)   
    616625                               ENDDO 
     
    623632                               k = surf_lsm_h%k(m)
    624633                               tend(k,j,i) = tend(k,j,i) + g / pt_reference    &
     634                                                         * drho_air_zw(k-1)    &
    625635                                                         * surf_lsm_h%shf(m)   
    626636                            ENDDO 
     
    632642                               k = surf_usm_h%k(m)
    633643                               tend(k,j,i) = tend(k,j,i) + g / pt_reference    &
     644                                                         * drho_air_zw(k-1)    &
    634645                                                         * surf_usm_h%shf(m)   
    635646                            ENDDO 
     
    641652                            DO  m = surf_s, surf_e
    642653                               k = surf_def_h(2)%k(m)
    643                                tend(k,j,i) = tend(k,j,i) + g / pt_reference *     &
    644                                                            surf_def_h(2)%shf(m)
     654                               tend(k,j,i) = tend(k,j,i) + g / pt_reference    &
     655                                                         * drho_air_zw(k-1)    &
     656                                                         * surf_def_h(2)%shf(m)
    645657                            ENDDO
    646658                         ENDIF
     
    754766                            tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *       &
    755767                                               ( k1 * surf_def_h(l)%shf(m) +   &
    756                                                  k2 * surf_def_h(l)%qsws(m) )
     768                                                 k2 * surf_def_h(l)%qsws(m)       &
     769                                               ) * drho_air_zw(k-1)
    757770                         ENDDO
    758771                      ENDDO
     
    788801                         tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *          &
    789802                                               ( k1 * surf_lsm_h%shf(m) +      &
    790                                                  k2 * surf_lsm_h%qsws(m) )
     803                                                 k2 * surf_lsm_h%qsws(m)       &
     804                                               ) * drho_air_zw(k-1)
    791805                      ENDDO
    792806!
     
    821835                         tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *          &
    822836                                               ( k1 * surf_usm_h%shf(m) +      &
    823                                                  k2 * surf_usm_h%qsws(m) )
     837                                                 k2 * surf_usm_h%qsws(m)       &
     838                                               ) * drho_air_zw(k-1)
    824839                      ENDDO
    825840
     
    861876                         tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *          &
    862877                                               ( k1 * surf_def_h(2)%shf(m) +   &
    863                                                  k2 * surf_def_h(2)%qsws(m) )
     878                                                 k2 * surf_def_h(2)%qsws(m)       &
     879                                               ) * drho_air_zw(k-1)
    864880
    865881                      ENDDO
     
    886902
    887903       USE arrays_3d,                                                          &
    888            ONLY:  ddzw, dd2zu, kh, km, prho, pt, q, ql, tend, u, v, vpt, w
     904           ONLY:  ddzw, dd2zu, drho_air_zw, kh, km, prho, pt, q, ql, tend, u,  &
     905                  v, vpt, w
    889906
    890907       USE cloud_parameters,                                                   &
     
    12541271                            k = surf_def_h(l)%k(m)
    12551272                            tend(k,j,i) = tend(k,j,i) + g / pt_reference *     &
    1256                                                            surf_def_h(l)%shf(m)
     1273                                                           drho_air_zw(k-1) *  &
     1274                                                           surf_def_h(l)%shf(m)
    12571275                         ENDDO
    12581276                      ENDDO
     
    12641282                         k = surf_lsm_h%k(m)
    12651283                         tend(k,j,i) = tend(k,j,i) + g / pt_reference *        &
     1284                                                           drho_air_zw(k-1) *  &
    12661285                                                           surf_lsm_h%shf(m)
    12671286                      ENDDO
     
    12731292                         k = surf_usm_h%k(m)
    12741293                         tend(k,j,i) = tend(k,j,i) + g / pt_reference *        &
     1294                                                           drho_air_zw(k-1) *  &
    12751295                                                           surf_usm_h%shf(m)
    12761296                      ENDDO
     
    12831303                         k = surf_def_h(2)%k(m)
    12841304                         tend(k,j,i) = tend(k,j,i) + g / pt_reference *        &
     1305                                                           drho_air_zw(k-1) *  &
    12851306                                                           surf_def_h(2)%shf(m)
    12861307                      ENDDO
     
    13311352                            k = surf_def_h(l)%k(m)
    13321353                            tend(k,j,i) = tend(k,j,i) + g / pt_reference       &
    1333                                                    * surf_def_h(l)%shf(m)   
     1354                                                   * drho_air_zw(k-1)          &
     1355                                                   * surf_def_h(l)%shf(m)
    13341356                         ENDDO 
    13351357                      ENDDO
     
    13411363                         k = surf_lsm_h%k(m)
    13421364                         tend(k,j,i) = tend(k,j,i) + g / pt_reference          &
    1343                                                      * surf_lsm_h%shf(m)   
     1365                                                     * drho_air_zw(k-1)        &
     1366                                                     * surf_lsm_h%shf(m)
    13441367                      ENDDO 
    13451368!
     
    13501373                         k = surf_usm_h%k(m)
    13511374                         tend(k,j,i) = tend(k,j,i) + g / pt_reference          &
    1352                                                      * surf_usm_h%shf(m)   
     1375                                                     * drho_air_zw(k-1)        &
     1376                                                     * surf_usm_h%shf(m)
    13531377                      ENDDO
    13541378                   ENDIF
     
    13601384                         k = surf_def_h(2)%k(m)
    13611385                         tend(k,j,i) = tend(k,j,i) + g / pt_reference *        &
    1362                                                            surf_def_h(2)%shf(m)
     1386                                                           drho_air_zw(k-1) *  &
     1387                                                           surf_def_h(2)%shf(m)
    13631388                      ENDDO
    13641389                   ENDIF
     
    14611486                      tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *             &
    14621487                                         ( k1 * surf_def_h(l)%shf(m) +         &
    1463                                            k2 * surf_def_h(l)%qsws(m) )
     1488                                           k2 * surf_def_h(l)%qsws(m)          &
     1489                                         ) * drho_air_zw(k-1)
    14641490                   ENDDO
    14651491                ENDDO
     
    14951521                   tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *                &
    14961522                                            ( k1 * surf_lsm_h%shf(m) +         &
    1497                                               k2 * surf_lsm_h%qsws(m) )
     1523                                              k2 * surf_lsm_h%qsws(m)          &
     1524                                            ) * drho_air_zw(k-1)
    14981525                ENDDO
    14991526!
     
    15281555                   tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *                &
    15291556                                            ( k1 * surf_usm_h%shf(m) +         &
    1530                                               k2 * surf_usm_h%qsws(m) )
     1557                                              k2 * surf_usm_h%qsws(m)          &
     1558                                            ) * drho_air_zw(k-1)
    15311559                ENDDO
    15321560
     
    15651593                   tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) *                &
    15661594                               ( k1* surf_def_h(2)%shf(m) +                    &
    1567                                  k2 * surf_def_h(2)%qsws(m) )
     1595                                 k2 * surf_def_h(2)%qsws(m)                    &
     1596                               ) * drho_air_zw(k-1)
    15681597                ENDDO
    15691598
     
    15851614
    15861615       USE arrays_3d,                                                          &
    1587            ONLY:  kh, km, u, v, zu
     1616           ONLY:  kh, km, drho_air_zw, u, v, zu
    15881617
    15891618       USE control_parameters,                                                 &
     
    16301659!--          effect of this error is negligible.
    16311660             surf_def_h(0)%u_0(m) = u(k+1,j,i) + surf_def_h(0)%usws(m) *       &
     1661                                        drho_air_zw(k-1) *                     &
    16321662                                        ( zu(k+1)    - zu(k-1)    )  /         &
    16331663                                        ( km(k,j,i)  + 1.0E-20_wp ) 
    16341664             surf_def_h(0)%v_0(m) = v(k+1,j,i) + surf_def_h(0)%vsws(m) *       &
     1665                                        drho_air_zw(k-1) *                     &
    16351666                                        ( zu(k+1)    - zu(k-1)    )  /         &
    16361667                                        ( km(k,j,i)  + 1.0E-20_wp ) 
     
    16611692!--          between u_0 and u(k-1).
    16621693             surf_def_h(1)%u_0(m) = u(k-1,j,i) - surf_def_h(1)%usws(m) *       &
     1694                                        drho_air_zw(k-1) *                     &
    16631695                                        ( zu(k+1)    - zu(k-1)    )  /         &
    16641696                                        ( km(k,j,i)  + 1.0E-20_wp ) 
    16651697             surf_def_h(1)%v_0(m) = v(k-1,j,i) - surf_def_h(1)%vsws(m) *       &
     1698                                        drho_air_zw(k-1) *                     &
    16661699                                        ( zu(k+1)    - zu(k-1)    )  /         &
    16671700                                        ( km(k,j,i)  + 1.0E-20_wp ) 
     
    16901723!--          effect of this error is negligible.
    16911724             surf_lsm_h%u_0(m) = u(k+1,j,i) + surf_lsm_h%usws(m)      *        &
     1725                                           drho_air_zw(k-1) *                  &
    16921726                                           ( zu(k+1)   - zu(k-1)    )  /       &
    16931727                                           ( km(k,j,i) + 1.0E-20_wp ) 
    16941728             surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m)      *        &
     1729                                           drho_air_zw(k-1) *                  &
    16951730                                           ( zu(k+1)   - zu(k-1)    )  /       &
    16961731                                           ( km(k,j,i) + 1.0E-20_wp )
     
    17191754!--          effect of this error is negligible.
    17201755             surf_usm_h%u_0(m) = u(k+1,j,i) + surf_usm_h%usws(m)      *        &
     1756                                           drho_air_zw(k-1) *                  &
    17211757                                           ( zu(k+1)   - zu(k-1)    )  /       &
    17221758                                           ( km(k,j,i) + 1.0E-20_wp ) 
    17231759             surf_usm_h%v_0(m) = v(k+1,j,i) + surf_usm_h%vsws(m)      *        &
     1760                                           drho_air_zw(k-1) *                  &
    17241761                                           ( zu(k+1)   - zu(k-1)    )  /       &
    17251762                                           ( km(k,j,i) + 1.0E-20_wp )
Note: See TracChangeset for help on using the changeset viewer.