Changeset 1007 for palm/trunk


Ignore:
Timestamp:
Sep 19, 2012 2:30:36 PM (12 years ago)
Author:
franke
Message:

Bugfixes


Missing calculation of mean particle weighting factor for output added. (data_output_2d, data_output_3d, data_output_mask, sum_up_3d_data)
Calculation of mean particle radius for output now considers the weighting factor. (data_output_mask)
Calculation of sugrid-scale buoyancy flux for humidity and cloud droplets corrected. (flow_statistics)
Factor in calculation of enhancement factor for collision efficencies corrected. (lpm_collision_kernels)
Calculation of buoyancy production now considers the liquid water mixing ratio in case of cloud droplets. (production_e)

Changes


Calculation of buoyancy flux for humidity in case of WS-scheme is now using turbulent fluxes of WS-scheme. (flow_statistics)
Calculation of the collision kernels now in SI units. (lpm_collision_kernels)

Location:
palm/trunk/SOURCE
Files:
7 edited

Legend:

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

    r979 r1007  
    44! Current revisions:
    55! -----------------
    6 !
     6! Bugfix: missing calculation of ql_vp added
    77!
    88! Former revisions:
     
    428428             CASE ( 'ql_vp_xy', 'ql_vp_xz', 'ql_vp_yz' )
    429429                IF ( av == 0 )  THEN
     430                   IF ( simulated_time >= particle_advection_start )  THEN
     431                      DO  i = nxl, nxr
     432                         DO  j = nys, nyn
     433                            DO  k = nzb, nzt+1
     434                               psi = prt_start_index(k,j,i)
     435                               DO  n = psi, psi+prt_count(k,j,i)-1
     436                                  tend(k,j,i) =  tend(k,j,i) + &
     437                                                 particles(n)%weight_factor / &
     438                                                 prt_count(k,j,i)
     439                               ENDDO
     440                            ENDDO
     441                         ENDDO
     442                      ENDDO
     443                      CALL exchange_horiz( tend, nbgp )
     444                   ELSE
     445                      tend = 0.0
     446                   END IF
     447                   DO  i = nxlg, nxrg
     448                      DO  j = nysg, nyng
     449                         DO  k = nzb, nzt+1
     450                            local_pf(i,j,k) = tend(k,j,i)
     451                         ENDDO
     452                      ENDDO
     453                   ENDDO
     454                   resorted = .TRUE.
     455                ELSE
     456                   CALL exchange_horiz( ql_vp_av, nbgp )
    430457                   to_be_resorted => ql_vp
    431                 ELSE
    432                    to_be_resorted => ql_vp_av
    433458                ENDIF
    434459                IF ( mode == 'xy' )  level_z = zu
  • palm/trunk/SOURCE/data_output_3d.f90

    r791 r1007  
    44! Current revisions:
    55! -----------------
    6 !
     6! Bugfix: missing calculation of ql_vp added
    77!
    88! Former revisions:
     
    295295          CASE ( 'ql_vp' )
    296296             IF ( av == 0 )  THEN
    297                 to_be_resorted => ql_vp
    298              ELSE
     297                DO  i = nxl, nxr
     298                   DO  j = nys, nyn
     299                      DO  k = nzb, nz_do3d
     300                         psi = prt_start_index(k,j,i)
     301                         DO  n = psi, psi+prt_count(k,j,i)-1
     302                            tend(k,j,i) = tend(k,j,i) + &
     303                                          particles(n)%weight_factor / &
     304                                          prt_count(k,j,i)
     305                         ENDDO
     306                      ENDDO
     307                   ENDDO
     308                ENDDO
     309                CALL exchange_horiz( tend, nbgp )
     310                DO  i = nxlg, nxrg
     311                   DO  j = nysg, nyng
     312                      DO  k = nzb, nz_do3d
     313                         local_pf(i,j,k) = tend(k,j,i)
     314                      ENDDO
     315                   ENDDO
     316                ENDDO
     317                resorted = .TRUE.
     318             ELSE
     319                CALL exchange_horiz( ql_vp_av, nbgp )
    299320                to_be_resorted => ql_vp_av
    300321             ENDIF
  • palm/trunk/SOURCE/data_output_mask.f90

    r772 r1007  
    44! Current revisions:
    55! -----------------
     6! Bugfix: calculation of pr must depend on the particle weighting factor,
     7! missing calculation of ql_vp added
    68!
    79! Former revisions:
     
    159161                         s_r4 = 0.0
    160162                         DO  n = psi, psi+prt_count(k,j,i)-1
    161                             s_r3 = s_r3 + particles(n)%radius**3
    162                             s_r4 = s_r4 + particles(n)%radius**4
     163                            s_r3 = s_r3 + particles(n)%radius**3 * &
     164                                          particles(n)%weight_factor
     165                            s_r4 = s_r4 + particles(n)%radius**4 * &
     166                                          particles(n)%weight_factor
    163167                         ENDDO
    164168                         IF ( s_r3 /= 0.0 )  THEN
     
    213217                to_be_resorted => q_av
    214218             ENDIF
    215              
     219
    216220          CASE ( 'ql' )
    217221             IF ( av == 0 )  THEN
     
    237241          CASE ( 'ql_vp' )
    238242             IF ( av == 0 )  THEN
    239                 to_be_resorted => ql_vp
    240              ELSE
     243                DO  i = nxl, nxr
     244                   DO  j = nys, nyn
     245                      DO  k = nzb, nz_do3d
     246                         psi = prt_start_index(k,j,i)
     247                         DO  n = psi, psi+prt_count(k,j,i)-1
     248                            tend(k,j,i) = tend(k,j,i) + &
     249                                          particles(n)%weight_factor / &
     250                                          prt_count(k,j,i)
     251                         ENDDO
     252                      ENDDO
     253                   ENDDO
     254                ENDDO
     255                CALL exchange_horiz( tend, nbgp )
     256                DO  i = 1, mask_size_l(mid,1)
     257                   DO  j = 1, mask_size_l(mid,2)
     258                      DO  k = 1, mask_size_l(mid,3)
     259                         local_pf(i,j,k) =  tend(mask_k(mid,k), &
     260                                   mask_j(mid,j),mask_i(mid,i))
     261                      ENDDO
     262                   ENDDO
     263                ENDDO
     264                resorted = .TRUE.
     265             ELSE
     266                CALL exchange_horiz( ql_vp_av, nbgp )
    241267                to_be_resorted => ql_vp_av
    242268             ENDIF
     
    264290                to_be_resorted => rho_av
    265291             ENDIF
    266              
     292
    267293          CASE ( 's' )
    268294             IF ( av == 0 )  THEN
     
    271297                to_be_resorted => s_av
    272298             ENDIF
    273              
     299
    274300          CASE ( 'sa' )
    275301             IF ( av == 0 )  THEN
     
    278304                to_be_resorted => sa_av
    279305             ENDIF
    280              
     306
    281307          CASE ( 'u' )
    282308             IF ( av == 0 )  THEN
  • palm/trunk/SOURCE/flow_statistics.f90

    r802 r1007  
    44! Current revisions:
    55! -----------------
     6! Calculation of buoyancy flux for humidity in case of WS-scheme is now using
     7! turbulent fluxes of WS-scheme
     8! Bugfix: Calculation of subgridscale buoyancy flux for humidity and cloud
     9! droplets at nzb and nzt added
    610!
    711! Former revisions:
     
    151155                        'timestep'
    152156       CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
    153      
     157
    154158    ENDIF
    155159
     
    176180       IF ( ws_scheme_mom .AND. sr == 0 )  THEN
    177181
    178 !       
     182!
    179183!--       According to the Neumann bc for the horizontal velocity components,
    180184!--       the corresponding fluxes has to satisfiy the same bc.
    181185          IF ( ocean )  THEN
    182186             sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
    183              sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)   
     187             sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)
    184188          ENDIF
    185189
    186190          DO  i = 0, threads_per_task-1
    187 !         
     191!
    188192!--          Swap the turbulent quantities evaluated in advec_ws.
    189193             sums_l(:,13,i) = sums_wsus_ws_l(:,i)       ! w*u*
     
    453457                                  ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)
    454458                ENDIF
    455                
     459
    456460!
    457461!--             Higher moments
     
    578582                                               * ( q(k+1,j,i) - q(k,j,i) )     &
    579583                                               * ddzu(k+1) * rmask(j,i,sr)
     584
    580585                   IF ( cloud_physics ) THEN
    581586                      sums_l(k,51,tn) = sums_l(k,51,tn)                        &
     
    618623                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
    619624                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
     625                   sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (           &
     626                                       ( 1.0 + 0.61 * q(nzb,j,i) ) *   &
     627                                       shf(j,i) + 0.61 * pt(nzb,j,i) * &
     628                                                  qsws(j,i) )
     629                   IF ( cloud_droplets )  THEN
     630                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (        &
     631                                         ( 1.0 + 0.61 * q(nzb,j,i) -   &
     632                                           ql(nzb,j,i) ) * shf(j,i) +  &
     633                                           0.61 * pt(nzb,j,i) * qsws(j,i) )
     634                   ENDIF
    620635                   IF ( cloud_physics )  THEN
    621                       sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (           &
    622                                           ( 1.0 + 0.61 * q(nzb,j,i) ) *   &
    623                                           shf(j,i) + 0.61 * pt(nzb,j,i) * &
    624                                                      qsws(j,i)            &
    625                                                               )
    626636!
    627637!--                   Formula does not work if ql(nzb) /= 0.0
     
    657667                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
    658668                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
     669                   sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (             &
     670                                       ( 1.0 + 0.61 * q(nzt,j,i) ) *     &
     671                                       tswst(j,i) + 0.61 * pt(nzt,j,i) * &
     672                                                  qswst(j,i) )
     673                   IF ( cloud_droplets )  THEN
     674                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (           &
     675                                          ( 1.0 + 0.61 * q(nzt,j,i) -     &
     676                                            ql(nzt,j,i) ) * tswst(j,i) +  &
     677                                            0.61 * pt(nzt,j,i) * qswst(j,i) )
     678                   ENDIF
    659679                   IF ( cloud_physics )  THEN
    660                       sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (           &
    661                                           ( 1.0 + 0.61 * q(nzt,j,i) ) *   &
    662                                           tswst(j,i) + 0.61 * pt(nzt,j,i) * &
    663                                                      qsws(j,i)            &
    664                                                               )
    665680!
    666681!--                   Formula does not work if ql(nzb) /= 0.0
     
    714729!--             content
    715730                IF ( humidity )  THEN
    716                    pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) + &
    717                                  vpt(k+1,j,i) - hom(k+1,1,44,sr) )
    718                    sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
    719                                                        rmask(j,i,sr)
    720 
    721                    IF ( cloud_physics  .OR.  cloud_droplets )  THEN
    722                       pts = 0.5 *                                            &
    723                              ( ( q(k,j,i)   - ql(k,j,i)   ) - hom(k,1,42,sr) &
    724                              + ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) )
    725                       sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) * &
     731                   IF ( cloud_physics .OR. cloud_droplets )  THEN
     732                      pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) + &
     733                                    vpt(k+1,j,i) - hom(k+1,1,44,sr) )
     734                      sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
    726735                                                          rmask(j,i,sr)
    727736                      sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * &
    728737                                                          rmask(j,i,sr)
    729                    ENDIF
    730                 ENDIF
    731 
     738                   ELSE
     739                      IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
     740                         pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) + &
     741                                       vpt(k+1,j,i) - hom(k+1,1,44,sr) )
     742                         sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
     743                                                             rmask(j,i,sr)
     744                      ELSE IF ( ws_scheme_sca .AND. sr == 0 )  THEN
     745                         sums_l(k,46,tn) = ( 1.0 + 0.61 * hom(k,1,41,sr) ) * &
     746                                             sums_l(k,17,tn) +      &
     747                                          0.61 * hom(k,1,4,sr) * sums_l(k,49,tn)
     748                      END IF
     749                   END IF
     750                ENDIF
    732751!
    733752!--             Passive scalar flux
     
    762781                  vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
    763782                              v(k+1,j,i) - hom(k+1,1,2,sr) )
    764 !                             
     783!
    765784!--               Momentum flux w*u*
    766785                  sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5 *                   &
     
    10021021          ENDDO
    10031022       ENDIF
     1023
    10041024!
    10051025!--    Collect horizontal average in hom.
  • palm/trunk/SOURCE/lpm_collision_kernels.f90

    r850 r1007  
    44! Current revisions:
    55! -----------------
    6 !
     6! converted all units to SI units and replaced some parameters by corresponding
     7! PALM parameters
     8! Bugfix: factor in calculation of enhancement factor for collision efficencies
     9! changed from 10. to 1.0
    710!
    811! Former revisions:
     
    6467
    6568    PUBLIC  ckernel, collision_efficiency_rogers, init_kernels, &
    66             rclass_lbound, rclass_ubound, recalculate_kernel           
     69            rclass_lbound, rclass_ubound, recalculate_kernel
    6770
    6871    REAL ::  epsilon, eps2, rclass_lbound, rclass_ubound, urms, urms2
     
    126129!             ENDIF
    127130          ENDDO
    128 !
    129 !--       Collision routines expect radius to be in cm
    130           radclass = radclass * 100.0
    131 
    132 !
    133 !--       Set the class bounds for dissipation in interval [0, 1000] cm**2/s**3
     131
     132!
     133!--       Set the class bounds for dissipation in interval [0.0, 0.1] m**2/s**3
    134134          DO  i = 1, dissipation_classes
    135              epsclass(i) = 1000.0 * REAL( i ) / dissipation_classes
     135             epsclass(i) = 0.1 * REAL( i ) / dissipation_classes
    136136!             IF ( myid == 0 )  THEN
    137137!                PRINT*, 'i=', i, ' eps = ', epsclass(i)
     
    148148
    149149             epsilon = epsclass(k)
    150              urms    = 202.0 * ( epsilon / 400.0 )**( 1.0 / 3.0 )
     150             urms    = 2.02 * ( epsilon / 0.04 )**( 1.0 / 3.0 )
    151151
    152152             CALL turbsd
     
    183183
    184184             PRINT*, '*** Hall kernel'
    185              WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E4, i = 1,radius_classes )
     185             WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E6, &
     186                                              i = 1,radius_classes )
    186187             DO  j = 1, radius_classes
    187                 WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j), ( hkernel(i,j), i = 1,radius_classes )
     188                WRITE ( *,'(F4.0,1X,20(F8.4,1X))' ) radclass(j),  &
     189                                          ( hkernel(i,j), i = 1,radius_classes )
    188190             ENDDO
    189191
     
    200202
    201203                PRINT*, '*** epsilon = ', epsclass(k)
    202                 WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E4, i = 1,radius_classes )
     204                WRITE ( *,'(5X,20(F4.0,1X))' ) ( radclass(i)*1.0E6, &
     205                                                 i = 1,radius_classes )
    203206                DO  j = 1, radius_classes
    204 !                   WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)*1.0E4, ( ckernel(i,j,k), i = 1,radius_classes )
    205                    WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)*1.0E4, ( hwratio(i,j), i = 1,radius_classes )
     207!                   WRITE ( *,'(F4.0,1X,20(F4.2,1X))' ) radclass(j)*1.0E6, &
     208!                                       ( ckernel(i,j,k), i = 1,radius_classes )
     209                   WRITE ( *,'(F4.0,1X,20(F8.4,1X))' ) radclass(j)*1.0E6, &
     210                                          ( hwratio(i,j), i = 1,radius_classes )
    206211                ENDDO
    207212             ENDDO
     
    210215
    211216          DEALLOCATE( ec, ecf, epsclass, gck, hkernel, winf )
    212 
    213           ckernel = ckernel * 1.0E-6   ! kernel is needed in m**3/s
    214217
    215218       ELSEIF( collision_kernel == 'hall'  .OR.  collision_kernel == 'wang' ) &
     
    249252
    250253!
    251 !--    Store particle radii on the radclass array. Collision routines
    252 !--    expect radii to be in cm.
    253        radclass(1:radius_classes) = particles(pstart:pend)%radius * 100.0
     254!--    Store particle radii on the radclass array
     255       radclass(1:radius_classes) = particles(pstart:pend)%radius
    254256
    255257       IF ( wang_kernel )  THEN
    256           epsilon = diss(k1,j1,i1) * 1.0E4   ! dissipation rate in cm**2/s**-3
     258          epsilon = diss(k1,j1,i1)   ! dissipation rate in m**2/s**3
    257259       ELSE
    258260          epsilon = 0.0
    259261       ENDIF
    260        urms    = 202.0 * ( epsilon / 400.0 )**( 0.33333333333 )
    261 
    262        IF ( wang_kernel  .AND.  epsilon > 0.001 )  THEN
     262       urms    = 2.02 * ( epsilon / 0.04 )**( 0.33333333333 )
     263
     264       IF ( wang_kernel  .AND.  epsilon > 1.0E-7 )  THEN
    263265!
    264266!--       Call routines to calculate efficiencies for the Wang kernel
     
    294296       ENDIF
    295297
    296        ckernel = ckernel * 1.0E-6   ! kernel is needed in m**3/s
    297 
    298298       DEALLOCATE( ec, radclass, winf )
    299299
     
    314314       USE particle_attributes
    315315       USE arrays_3d
     316       USE control_parameters
    316317
    317318       IMPLICIT NONE
     
    326327                v1xysq, v2, v2xysq, wrfin, wrgrav2, wrtur2xy, xx, yy, z
    327328
    328        REAL, SAVE ::  airdens, airvisc, anu, gravity, waterdens
    329 
    330329       REAL, DIMENSION(1:radius_classes) ::  st, tau
    331330
     
    336335
    337336          first = .FALSE.
    338           airvisc   = 0.1818     ! dynamic viscosity in mg/cm*s
    339           airdens   = 1.2250     ! air density in mg/cm**3
    340           waterdens = 1000.0     ! water density in mg/cm**3
    341           gravity   = 980.6650   ! in cm/s**2
    342           anu = airvisc/airdens  ! kinetic viscosity in cm**2/s
    343 
    344        ENDIF   
    345 
    346        lambda    = urms * SQRT( 15.0 * anu / epsilon )     ! in cm
    347        lambda_re = urms**2 * SQRT( 15.0 / epsilon / anu )
     337
     338       ENDIF
     339
     340       lambda    = urms * SQRT( 15.0 * molecular_viscosity / epsilon )    ! in m
     341       lambda_re = urms**2 * SQRT( 15.0 / epsilon / molecular_viscosity )
    348342       tl        = urms**2 / epsilon                       ! in s
    349        lf        = 0.5 * urms**3 / epsilon                 ! in cm
    350        tauk      = SQRT( anu / epsilon )                   ! in s
    351        eta       = ( anu**3 / epsilon )**0.25              ! in cm
    352        vk        = eta / tauk                              ! in cm/s
     343       lf        = 0.5 * urms**3 / epsilon                 ! in m
     344       tauk      = SQRT( molecular_viscosity / epsilon )                  ! in s
     345       eta       = ( molecular_viscosity**3 / epsilon )**0.25             ! in m
     346       vk        = eta / tauk
    353347
    354348       ao = ( 11.0 + 7.0 * lambda_re ) / ( 205.0 + lambda_re )
    355349       tt = SQRT( 2.0 * lambda_re / ( SQRT( 15.0 ) * ao ) ) * tauk   ! in s
    356350
    357        CALL fallg    ! gives winf in cm/s
     351       CALL fallg    ! gives winf in m/s
    358352
    359353       DO  i = 1, radius_classes
    360           tau(i) = winf(i) / gravity    ! in s
     354          tau(i) = winf(i) / g    ! in s
    361355          st(i)  = tau(i) / tauk
    362356       ENDDO
     
    368362       bbb = SQRT( 1.0 - 2.0 * be**2 )
    369363       d1  = ( 1.0 + bbb ) / ( 2.0 * bbb )
    370        e1  = lf * ( 1.0 + bbb ) * 0.5   ! in cm
     364       e1  = lf * ( 1.0 + bbb ) * 0.5   ! in m
    371365       d2  = ( 1.0 - bbb ) * 0.5 / bbb
    372        e2  = lf * ( 1.0 - bbb ) * 0.5   ! in cm
     366       e2  = lf * ( 1.0 - bbb ) * 0.5   ! in m
    373367       ccc = SQRT( 1.0 - 2.0 * z**2 )
    374368       b1  = ( 1.0 + ccc ) * 0.5 / ccc
     
    379373       DO  i = 1, radius_classes
    380374
    381           v1 = winf(i)        ! in cm/s
     375          v1 = winf(i)        ! in m/s
    382376          t1 = tau(i)         ! in s
    383377
    384378          DO  j = 1, i
    385              rrp = radclass(i) + radclass(j)               ! radius in cm
    386              v2  = winf(j)                                 ! in cm/s
     379             rrp = radclass(i) + radclass(j)
     380             v2  = winf(j)                                 ! in m/s
    387381             t2  = tau(j)                                  ! in s
    388382
    389              v1xysq  = b1 * d1 * phi(c1,e1,v1,t1) - b1 * d2 * phi(c1,e2,v1,t1) &
    390                      - b2 * d1 * phi(c2,e1,v1,t1) + b2 * d2 * phi(c2,e2,v1,t1)
    391              v1xysq  = v1xysq * urms**2 / t1                ! in cm**2/s**2
    392              vrms1xy = SQRT( v1xysq )                       ! in cm/s
    393 
    394              v2xysq  = b1 * d1 * phi(c1,e1,v2,t2) - b1 * d2 * phi(c1,e2,v2,t2) &
    395                      - b2 * d1 * phi(c2,e1,v2,t2) + b2 * d2 * phi(c2,e2,v2,t2)
    396              v2xysq  = v2xysq * urms**2 / t2                ! in cm**2/s**2
    397              vrms2xy = SQRT( v2xysq )                       ! in cm/s
     383             v1xysq  = b1 * d1 * phi_w(c1,e1,v1,t1) - b1 * d2 * phi_w(c1,e2,v1,t1) &
     384                     - b2 * d1 * phi_w(c2,e1,v1,t1) + b2 * d2 * phi_w(c2,e2,v1,t1)
     385             v1xysq  = v1xysq * urms**2 / t1                ! in m**2/s**2
     386             vrms1xy = SQRT( v1xysq )                       ! in m/s
     387
     388             v2xysq  = b1 * d1 * phi_w(c1,e1,v2,t2) - b1 * d2 * phi_w(c1,e2,v2,t2) &
     389                     - b2 * d1 * phi_w(c2,e1,v2,t2) + b2 * d2 * phi_w(c2,e2,v2,t2)
     390             v2xysq  = v2xysq * urms**2 / t2                ! in m**2/s**2
     391             vrms2xy = SQRT( v2xysq )                       ! in m/s
    398392
    399393             IF ( winf(i) >= winf(j) )  THEN
     
    414408                         b2 * d2* zhi(c2,e2,v1,t1,v2,t2)
    415409             fr       = d1 * EXP( -rrp / e1 ) - d2 * EXP( -rrp / e2 )
    416              v1v2xy   = v1v2xy * fr * urms**2 / tau(i) / tau(j)  ! in cm**2/s**2
    417              wrtur2xy = vrms1xy**2 + vrms2xy**2 - 2.0 * v1v2xy   ! in cm**2/s**2
     410             v1v2xy   = v1v2xy * fr * urms**2 / tau(i) / tau(j)  ! in m**2/s**2
     411             wrtur2xy = vrms1xy**2 + vrms2xy**2 - 2.0 * v1v2xy   ! in m**2/s**2
    418412             IF ( wrtur2xy < 0.0 )  wrtur2xy = 0.0
    419413             wrgrav2  = pi / 8.0 * ( winf(j) - winf(i) )**2
    420              wrfin    = SQRT( ( 2.0 / pi ) * ( wrtur2xy + wrgrav2) )   ! in cm/s
     414             wrfin    = SQRT( ( 2.0 / pi ) * ( wrtur2xy + wrgrav2) )   ! in m/s
    421415
    422416!
     
    433427             yy = 0.1886 * EXP( 20.306 / lambda_re )
    434428
    435              c1_gr  =  xx / ( gravity / vk * tauk )**yy
    436 
    437              ao_gr  = ao + ( pi / 8.0) * ( gravity / vk * tauk )**2
     429             c1_gr  =  xx / ( g / vk * tauk )**yy
     430
     431             ao_gr  = ao + ( pi / 8.0) * ( g / vk * tauk )**2
    438432             fao_gr = 20.115 * SQRT( ao_gr / lambda_re )
    439433             rc     = SQRT( fao_gr * ABS( st(j) - st(i) ) ) * eta   ! in cm
     
    452446
    453447!------------------------------------------------------------------------------!
    454 ! phi as a function
    455 !------------------------------------------------------------------------------!
    456     REAL FUNCTION phi( a, b, vsett, tau0 )
     448! phi_w as a function
     449!------------------------------------------------------------------------------!
     450    REAL FUNCTION phi_w( a, b, vsett, tau0 )
    457451
    458452       IMPLICIT NONE
     
    461455
    462456       aa1 = 1.0 / tau0 + 1.0 / a + vsett / b
    463        phi = 1.0 / aa1  - 0.5 * vsett / b / aa1**2  ! in s
    464 
    465     END FUNCTION phi
    466 
    467 
    468 !------------------------------------------------------------------------------!
    469 ! zeta as a function
     457       phi_w = 1.0 / aa1  - 0.5 * vsett / b / aa1**2  ! in s
     458
     459    END FUNCTION phi_w
     460
     461
     462!------------------------------------------------------------------------------!
     463! zhi as a function
    470464!------------------------------------------------------------------------------!
    471465    REAL FUNCTION zhi( a, b, vsett1, tau1, vsett2, tau2 )
     
    490484
    491485!------------------------------------------------------------------------------!
    492 ! Calculation of terminal velocity winf
     486! Calculation of terminal velocity winf following Equations 10-138 to 10-145
     487! from (Pruppacher and Klett, 1997)
    493488!------------------------------------------------------------------------------!
    494489    SUBROUTINE fallg
     
    498493       USE particle_attributes
    499494       USE arrays_3d
     495       USE control_parameters
    500496
    501497       IMPLICIT NONE
     
    505501       LOGICAL, SAVE ::  first = .TRUE.
    506502
    507        REAL, SAVE ::  cunh, eta, grav, phy, py, rhoa, rhow, sigma, stb, stok, &
     503       REAL, SAVE ::  cunh, eta, phy, py, rho_a, sigma, stb, stok, &
    508504                      t0, xlamb
    509505
     
    523519                 -0.542819E-1, 0.238449E-2 /)
    524520
    525           eta   = 1.818E-4          ! in poise = g/(cm s)
    526           xlamb = 6.62E-6           ! in cm
    527           rhow  = 1.0               ! in g/cm**3
    528           rhoa  = 1.225E-3          ! in g/cm**3
    529           grav  = 980.665           ! in cm/s**2
    530           cunh  = 1.257 * xlamb     ! in cm
    531           t0    = 273.15            ! in K
    532           sigma = 76.1 - 0.155 * ( 293.15 - t0 )              ! in N/m = g/s**2
    533           stok  = 2.0  * grav * ( rhow - rhoa ) / ( 9.0 * eta ) ! in 1/(cm s)
    534           stb   = 32.0 * rhoa * ( rhow - rhoa) * grav / (3.0 * eta * eta)
    535                                                                 ! in 1/cm**3
    536           phy   = sigma**3 * rhoa**2 / ( eta**4 * grav * ( rhow - rhoa ) )
     521!
     522!--       Parameter values for p = 1013,25 hPa and T = 293,15 K
     523          eta   = 1.818E-5         ! in kg/(m s)
     524          xlamb = 6.6E-8           ! in m
     525          rho_a = 1.204            ! in kg/m**3
     526          cunh  = 1.26 * xlamb     ! in m
     527          sigma = 0.07363          ! in kg/s**2
     528          stok  = 2.0  * g * ( rho_l - rho_a ) / ( 9.0 * eta ) ! in 1/(m s)
     529          stb   = 32.0 * rho_a * ( rho_l - rho_a) * g / (3.0 * eta * eta)
     530          phy   = sigma**3 * rho_a**2 / ( eta**4 * g * ( rho_l - rho_a ) )
    537531          py    = phy**( 1.0 / 6.0 )
    538532
     
    541535       DO  j = 1, radius_classes
    542536
    543           IF ( radclass(j) <= 1.0E-3 ) THEN
    544 
    545              winf(j) = stok * ( radclass(j)**2 + cunh * radclass(j) ) ! in cm/s
    546 
    547           ELSEIF ( radclass(j) > 1.0E-3  .AND.  radclass(j) <= 5.35E-2 )  THEN
     537          IF ( radclass(j) <= 1.0E-5 ) THEN
     538
     539             winf(j) = stok * ( radclass(j)**2 + cunh * radclass(j) )
     540
     541          ELSEIF ( radclass(j) > 1.0E-5  .AND.  radclass(j) <= 5.35E-4 )  THEN
    548542
    549543             x = LOG( stb * radclass(j)**3 )
     
    553547                y = y + b(i) * x**(i-1)
    554548             ENDDO
    555 
    556              xrey = ( 1.0 + cunh / radclass(j) ) * EXP( y )
    557              winf(j) = xrey * eta / ( 2.0 * rhoa * radclass(j) )      ! in cm/s
    558 
    559           ELSEIF ( radclass(j) > 5.35E-2 )  THEN
    560 
    561              IF ( radclass(j) > 0.35 )  THEN
    562                 bond = grav * ( rhow - rhoa ) * 0.35**2 / sigma
     549!
     550!--          Note: this Eq. is wrong in (Pruppacher and Klett, 1997, p. 418)
     551!--          for correct version see (Beard, 1976)
     552             xrey = ( 1.0 + cunh / radclass(j) ) * EXP( y )
     553
     554             winf(j) = xrey * eta / ( 2.0 * rho_a * radclass(j) )
     555
     556          ELSEIF ( radclass(j) > 5.35E-4 )  THEN
     557
     558             IF ( radclass(j) > 0.0035 )  THEN
     559                bond = g * ( rho_l - rho_a ) * 0.0035**2 / sigma
    563560             ELSE
    564                 bond = grav * ( rhow - rhoa ) * radclass(j)**2 / sigma
     561               bond = g * ( rho_l - rho_a ) * radclass(j)**2 / sigma
    565562             ENDIF
    566563
     
    574571             xrey = py * EXP( y )
    575572
    576              IF ( radclass(j) > 0.35 )  THEN
    577                 winf(j) = xrey * eta / ( 2.0 * rhoa * 0.35 )           ! in cm/s
     573             IF ( radclass(j) > 0.0035 )  THEN
     574                winf(j) = xrey * eta / ( 2.0 * rho_a * 0.0035 )
    578575             ELSE
    579                 winf(j) = xrey * eta / ( 2.0 * rhoa * radclass(j) )    ! in cm/s
     576                winf(j) = xrey * eta / ( 2.0 * rho_a * radclass(j) )
    580577             ENDIF
    581578
     
    668665!
    669666!--    Calculate the radius class index of particles with respect to array r
     667!--    Radius has to be in µm
    670668       ALLOCATE( ira(1:radius_classes) )
    671669       DO  j = 1, radius_classes
    672           particle_radius = radclass(j) * 1.0E4   ! in microm
     670          particle_radius = radclass(j) * 1.0E6
    673671          DO  k = 1, 15
    674672             IF ( particle_radius < r0(k) )  THEN
     
    693691             IF ( ir < 16 )  THEN
    694692                IF ( ir >= 2 )  THEN
    695                    pp = ( ( radclass(j) * 1.0E04 ) - r0(ir-1) ) / &
     693                   pp = ( ( radclass(j) * 1.0E06 ) - r0(ir-1) ) / &
    696694                        ( r0(ir) - r0(ir-1) )
    697695                   qq = ( rq- rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
     
    754752          rat = (/ 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0 /)
    755753!
    756 !--       In 100 cm^2/s^3
     754!--       for 100 cm**2/s**3
    757755          ecoll_100(:,1) = (/1.74,  1.74,  1.773, 1.49,  1.207,  1.207,  1.0 /)
    758756          ecoll_100(:,2) = (/1.46,  1.46,  1.421, 1.245, 1.069,  1.069,  1.0 /)
     
    767765          ecoll_100(:,11)= (/20.3,  20.3,  14.6 , 8.61,  2.60,   2.60 ,  1.0 /)
    768766!
    769 !--       In 400 cm^2/s^3
     767!--       for 400 cm**2/s**3
    770768          ecoll_400(:,1) = (/4.976, 4.976,  3.593, 2.519, 1.445,  1.445,  1.0 /)
    771769          ecoll_400(:,2) = (/2.984, 2.984,  2.181, 1.691, 1.201,  1.201,  1.0 /)
     
    784782!
    785783!--    Calculate the radius class index of particles with respect to array r0
     784!--    Radius has to be in µm
    786785       ALLOCATE( ira(1:radius_classes) )
    787786
    788787       DO  j = 1, radius_classes
    789           particle_radius = radclass(j) * 1.0E4    ! in microm
     788          particle_radius = radclass(j) * 1.0E6
    790789          DO  k = 1, 7
    791790             IF ( particle_radius < r0(k) )  THEN
     
    799798!
    800799!--    Two-dimensional linear interpolation of the collision efficiencies
     800!--    Radius has to be in µm
    801801       DO  j =  1, radius_classes
    802802          DO  i = 1, j
     
    812812             ENDDO
    813813
    814              y1 = 1.0      ! 0  cm2/s3
    815 !
    816 !--          100 cm2/s3, 400 cm2/s3
     814             y1 = 0.0001      ! for 0 m**2/s**3
     815
    817816             IF ( ir < 8 )  THEN
    818817                IF ( ir >= 2 )  THEN
    819                    pp = ( radclass(j)*1.0E4 - r0(ir-1) ) / ( r0(ir) - r0(ir-1) )
     818                   pp = ( radclass(j)*1.0E6 - r0(ir-1) ) / ( r0(ir) - r0(ir-1) )
    820819                   qq = ( rq - rat(iq-1) ) / ( rat(iq) - rat(iq-1) )
    821820                   y2 = ( 1.0-pp ) * ( 1.0-qq ) * ecoll_100(ir-1,iq-1) +  &
    822821                                pp * ( 1.0-qq ) * ecoll_100(ir,iq-1)   +  &
    823                                 qq * ( 10.-pp ) * ecoll_100(ir-1,iq)   +  &
     822                                qq * ( 1.0-pp ) * ecoll_100(ir-1,iq)   +  &
    824823                                pp * qq         * ecoll_100(ir,iq)
    825824                   y3 = ( 1.0-pp ) * ( 1.0-qq ) * ecoll_400(ir-1,iq-1) +  &
     
    838837             ENDIF
    839838!
    840 !--          Linear interpolation of dissipation rate in cm2/s3
    841              IF ( epsilon <= 100.0 )  THEN
    842                 ecf(j,i) = ( epsilon - 100.0 ) / (   0.0 - 100.0 ) * y1 &
    843                          + ( epsilon -   0.0 ) / ( 100.0 -   0.0 ) * y2
    844              ELSEIF ( epsilon <= 600.0 )  THEN
    845                 ecf(j,i) = ( epsilon - 400.0 ) / ( 100.0 - 400.0 ) * y2 &
    846                          + ( epsilon - 100.0 ) / ( 400.0 - 100.0 ) * y3
     839!--          Linear interpolation of dissipation rate in m**2/s**3
     840             IF ( epsilon <= 0.01 )  THEN
     841                ecf(j,i) = ( epsilon - 0.01 ) / (   0.0 - 0.01 ) * y1 &
     842                         + ( epsilon -   0.0 ) / ( 0.01 -   0.0 ) * y2
     843             ELSEIF ( epsilon <= 0.06 )  THEN
     844                ecf(j,i) = ( epsilon - 0.04 ) / ( 0.01 - 0.04 ) * y2 &
     845                         + ( epsilon - 0.01 ) / ( 0.04 - 0.01 ) * y3
    847846             ELSE
    848                 ecf(j,i) = (   600.0 - 400.0 ) / ( 100.0 - 400.0 ) * y2 &
    849                          + (   600.0 - 100.0 ) / ( 400.0 - 100.0 ) * y3
     847                ecf(j,i) = (   0.06 - 0.04 ) / ( 0.01 - 0.04 ) * y2 &
     848                         + (   0.06 - 0.01 ) / ( 0.04 - 0.01 ) * y3
    850849             ENDIF
    851850
  • palm/trunk/SOURCE/production_e.f90

    r941 r1007  
    44! Current revisions:
    55! -----------------
    6 !
     6! Bugfix: calculation of buoyancy production has to consider the liquid water
     7! mixing ratio in case of cloud droplets
    78!
    89! Former revisions:
     
    164165
    165166                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
    166                
     167
    167168             ENDDO
    168169          ENDDO
     
    189190
    190191                   IF ( wall_e_y(j,i) /= 0.0 )  THEN
    191 !                     
     192!
    192193!--                   Inconsistency removed: as the thermal stratification is
    193194!--                   not taken into account for the evaluation of the wall
     
    219220
    220221                   IF ( wall_e_x(j,i) /= 0.0 )  THEN
    221 !                     
     222!
    222223!--                   Inconsistency removed: as the thermal stratification is
    223224!--                   not taken into account for the evaluation of the wall
     
    275276
    276277                      IF ( wall_e_y(j,i) /= 0.0 )  THEN
    277 !                     
     278!
    278279!--                      Inconsistency removed: as the thermal stratification
    279280!--                      is not taken into account for the evaluation of the
     
    301302
    302303                      IF ( wall_e_x(j,i) /= 0.0 )  THEN
    303 !                     
     304!
    304305!--                      Inconsistency removed: as the thermal stratification
    305306!--                      is not taken into account for the evaluation of the
     
    416417
    417418                   tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
    418                
     419
    419420                ENDIF
    420421
     
    551552                   DO  k = nzb_diff_s_inner(j,i), nzt_diff
    552553
    553                       IF ( .NOT. cloud_physics ) THEN
     554                      IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
    554555                         k1 = 1.0 + 0.61 * q(k,j,i)
    555556                         k2 = 0.61 * pt(k,j,i)
    556                       ELSE
     557                         tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
     558                                         g / vpt(k,j,i) *                      &
     559                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
     560                                           k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
     561                                         ) * dd2zu(k)
     562                      ELSE IF ( cloud_physics )  THEN
    557563                         IF ( ql(k,j,i) == 0.0 )  THEN
    558564                            k1 = 1.0 + 0.61 * q(k,j,i)
     
    568574                            k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
    569575                         ENDIF
    570                       ENDIF
    571 
    572                       tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
     576                         tend(k,j,i) = tend(k,j,i) - kh(k,j,i) *               &
     577                                         g / vpt(k,j,i) *                      &
    573578                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +  &
    574579                                           k2 * ( q(k+1,j,i) - q(k-1,j,i) )    &
    575580                                         ) * dd2zu(k)
     581                      ELSE IF ( cloud_droplets )  THEN
     582                         k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
     583                         k2 = 0.61 * pt(k,j,i)
     584                         tend(k,j,i) = tend(k,j,i) -                          &
     585                                       kh(k,j,i) * g / vpt(k,j,i) *           &
     586                                       ( k1 * ( pt(k+1,j,i)- pt(k-1,j,i) ) +  &
     587                                         k2 * ( q(k+1,j,i) -  q(k-1,j,i) ) -  &
     588                                         pt(k,j,i) * ( ql(k+1,j,i) -          &
     589                                         ql(k-1,j,i) ) ) * dd2zu(k)
     590                      ENDIF
     591
    576592                   ENDDO
    577593
     
    584600                      k = nzb_diff_s_inner(j,i)-1
    585601
    586                       IF ( .NOT. cloud_physics ) THEN
     602                      IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
    587603                         k1 = 1.0 + 0.61 * q(k,j,i)
    588604                         k2 = 0.61 * pt(k,j,i)
    589                       ELSE
     605                      ELSE IF ( cloud_physics )  THEN
    590606                         IF ( ql(k,j,i) == 0.0 )  THEN
    591607                            k1 = 1.0 + 0.61 * q(k,j,i)
     
    601617                            k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
    602618                         ENDIF
     619                      ELSE IF ( cloud_droplets )  THEN
     620                         k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
     621                         k2 = 0.61 * pt(k,j,i)
    603622                      ENDIF
    604623
     
    615634                      k = nzt
    616635
    617                       IF ( .NOT. cloud_physics ) THEN
     636                      IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
    618637                         k1 = 1.0 + 0.61 * q(k,j,i)
    619638                         k2 = 0.61 * pt(k,j,i)
    620                       ELSE
     639                      ELSE IF ( cloud_physics )  THEN
    621640                         IF ( ql(k,j,i) == 0.0 )  THEN
    622641                            k1 = 1.0 + 0.61 * q(k,j,i)
     
    632651                            k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
    633652                         ENDIF
     653                      ELSE IF ( cloud_droplets )  THEN
     654                         k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
     655                         k2 = 0.61 * pt(k,j,i)
    634656                      ENDIF
    635657
     
    699721
    700722          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
    701                
     723
    702724       ENDDO
    703725
     
    721743
    722744             IF ( wall_e_y(j,i) /= 0.0 )  THEN
    723 !                     
     745!
    724746!--             Inconsistency removed: as the thermal stratification
    725747!--             is not taken into account for the evaluation of the
     
    751773
    752774             IF ( wall_e_x(j,i) /= 0.0 )  THEN
    753 !                     
     775!
    754776!--             Inconsistency removed: as the thermal stratification
    755777!--             is not taken into account for the evaluation of the
     
    805827
    806828                IF ( wall_e_y(j,i) /= 0.0 )  THEN
    807 !                     
     829!
    808830!--                Inconsistency removed: as the thermal stratification
    809831!--                is not taken into account for the evaluation of the
     
    831853
    832854                IF ( wall_e_x(j,i) /= 0.0 )  THEN
    833 !                     
     855!
    834856!--                Inconsistency removed: as the thermal stratification
    835857!--                is not taken into account for the evaluation of the
     
    10411063             DO  k = nzb_diff_s_inner(j,i), nzt_diff
    10421064
    1043                 IF ( .NOT. cloud_physics )  THEN
     1065                IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets )  THEN
    10441066                   k1 = 1.0 + 0.61 * q(k,j,i)
    10451067                   k2 = 0.61 * pt(k,j,i)
    1046                 ELSE
     1068                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *   &
     1069                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
     1070                                           k2 * ( q(k+1,j,i) - q(k-1,j,i) )   &
     1071                                         ) * dd2zu(k)
     1072                ELSE IF ( cloud_physics )  THEN
    10471073                   IF ( ql(k,j,i) == 0.0 )  THEN
    10481074                      k1 = 1.0 + 0.61 * q(k,j,i)
     
    10581084                      k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
    10591085                   ENDIF
    1060                 ENDIF
    1061 
    1062                 tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *      &
     1086                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *   &
    10631087                                         ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
    10641088                                           k2 * ( q(k+1,j,i) - q(k-1,j,i) )   &
    10651089                                         ) * dd2zu(k)
     1090                ELSE IF ( cloud_droplets )  THEN
     1091                   k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
     1092                   k2 = 0.61 * pt(k,j,i)
     1093                   tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) *  &
     1094                                     ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) +    &
     1095                                       k2 * ( q(k+1,j,i) - q(k-1,j,i) ) -    &
     1096                                       pt(k,j,i) * ( ql(k+1,j,i) -           &
     1097                                                     ql(k-1,j,i) ) ) * dd2zu(k)
     1098                ENDIF
    10661099             ENDDO
    10671100
     
    10691102                k = nzb_diff_s_inner(j,i)-1
    10701103
    1071                 IF ( .NOT. cloud_physics ) THEN
     1104                IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
    10721105                   k1 = 1.0 + 0.61 * q(k,j,i)
    10731106                   k2 = 0.61 * pt(k,j,i)
    1074                 ELSE
     1107                ELSE IF ( cloud_physics )  THEN
    10751108                   IF ( ql(k,j,i) == 0.0 )  THEN
    10761109                      k1 = 1.0 + 0.61 * q(k,j,i)
     
    10861119                      k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
    10871120                   ENDIF
     1121                ELSE IF ( cloud_droplets )  THEN
     1122                   k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
     1123                   k2 = 0.61 * pt(k,j,i)
    10881124                ENDIF
    10891125
     
    10951131                k = nzt
    10961132
    1097                 IF ( .NOT. cloud_physics ) THEN
     1133                IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
    10981134                   k1 = 1.0 + 0.61 * q(k,j,i)
    10991135                   k2 = 0.61 * pt(k,j,i)
    1100                 ELSE
     1136                ELSE IF ( cloud_physics )  THEN
    11011137                   IF ( ql(k,j,i) == 0.0 )  THEN
    11021138                      k1 = 1.0 + 0.61 * q(k,j,i)
     
    11121148                      k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
    11131149                   ENDIF
     1150                ELSE IF ( cloud_droplets )  THEN
     1151                   k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
     1152                   k2 = 0.61 * pt(k,j,i)
    11141153                ENDIF
    11151154
  • palm/trunk/SOURCE/sum_up_3d_data.f90

    r979 r1007  
    44! Current revisions:
    55! -----------------
    6 !
     6! Bugfix in calculation of ql_vp
    77!
    88! Former revisions:
     
    388388
    389389          CASE ( 'ql_vp' )
    390              DO  i = nxlg, nxrg
    391                 DO  j = nysg, nyng
    392                    DO  k = nzb, nzt+1
    393                       ql_vp_av(k,j,i) = ql_vp_av(k,j,i) + ql_vp(k,j,i)
     390             DO  i = nxl, nxr
     391                DO  j = nys, nyn
     392                   DO  k = nzb, nzt+1
     393                      psi = prt_start_index(k,j,i)
     394                      DO  n = psi, psi+prt_count(k,j,i)-1
     395                         ql_vp_av(k,j,i) = ql_vp_av(k,j,i) + &
     396                                           particles(n)%weight_factor / &
     397                                           prt_count(k,j,i)
     398                      ENDDO
    394399                   ENDDO
    395400                ENDDO
Note: See TracChangeset for help on using the changeset viewer.