Changeset 1065 for palm


Ignore:
Timestamp:
Nov 22, 2012 5:42:36 PM (12 years ago)
Author:
hoffmann
Message:

cloud physics: rain sedimentation and turbulence effects

Location:
palm/trunk/SOURCE
Files:
11 edited

Legend:

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

    r1061 r1065  
    2020! Current revisions:
    2121! -----------------
     22! Bugfix: It is not allowed to use cloud_scheme = seifert_beheng without
     23!         precipitation in order to save computational resources.
    2224!
    2325!
     
    886888                        'loop_optimization = cache'
    887889       CALL message( 'check_parameters', 'PA0362', 1, 2, 0, 6, 0 )
     890    ENDIF
     891
     892    IF ( cloud_physics  .AND. icloud_scheme == 0  .AND.                      &
     893         .NOT. precipitation ) THEN
     894       message_string = 'cloud_scheme = seifert_beheng requires ' // &
     895                        'precipitation = .TRUE.'
     896       CALL message( 'check_parameters', 'PA0363', 1, 2, 0, 6, 0 )
    888897    ENDIF
    889898
  • palm/trunk/SOURCE/data_output_2d.f90

    r1054 r1065  
    2020! Current revisions:
    2121! -----------------
     22! Bugfix: Output of cross sections of ql
    2223!
    2324! Former revisions:
     
    488489                      DO  i = nxlg, nxrg
    489490                         DO  j = nysg, nyng
    490                             DO  k = nzb, nz_do3d
     491                            DO  k = nzb, nzt+1
    491492                               local_pf(i,j,k) = ql(k,j,i) + qr(k,j,i)
    492493                            ENDDO
     
    501502                      DO  i = nxlg, nxrg
    502503                         DO  j = nysg, nyng
    503                             DO  k = nzb, nz_do3d
     504                            DO  k = nzb, nzt+1
    504505                               local_pf(i,j,k) = ql_av(k,j,i) + qr_av(k,j,i)
    505506                            ENDDO
  • palm/trunk/SOURCE/diffusion_e.f90

    r1037 r1065  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Enabled the claculation of diss in case of turbulence = .TRUE. (parameterized
     23! effects of turbulence on autoconversion and accretion in two-moments cloud
     24! physics scheme).
    2325!
    2426! Former revisions:
     
    183185!--          Store dissipation if needed for calculating the sgs particle
    184186!--          velocities
    185              IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
     187             IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.               &
     188                  turbulence )  THEN
    186189                DO  j = nys, nyn
    187190                   DO  k = nzb_s_inner(j,i)+1, nzt
     
    254257!--          Store dissipation if needed for calculating the sgs particle
    255258!--          velocities
    256              IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
     259             IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.               &
     260                  turbulence )  THEN
    257261                DO  j = nys, nyn
    258262                   DO  k = nzb_s_inner(j,i)+1, nzt
     
    268272!
    269273!--    Boundary condition for dissipation
    270        IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
     274       IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.  turbulence )  THEN
    271275          DO  i = nxl, nxr
    272276             DO  j = nys, nyn
     
    433437!--                   Store dissipation if needed for calculating the sgs
    434438!--                   particle  velocities
    435                       IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
     439                      IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.      &
     440                           turbulence )  THEN
    436441                         diss(k,j,i) = dissipation
    437442                      ENDIF
     
    448453!
    449454!--    Boundary condition for dissipation
    450        IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
     455       IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.  turbulence )  THEN
    451456          !$acc kernels present( diss, nzb_s_inner )
    452457          !$acc loop
     
    539544!
    540545!--    Store dissipation if needed for calculating the sgs particle velocities
    541        IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
     546       IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.  turbulence )  THEN
    542547          DO  k = nzb_s_inner(j,i)+1, nzt
    543548             diss(k,j,i) = dissipation(k)
  • palm/trunk/SOURCE/init_3d_model.f90

    r1054 r1065  
    2323! Current revisions:
    2424! ------------------
     25! allocation of diss (dissipation rate) in case of turbulence = .TRUE. added
    2526!
    2627! Former revisions:
     
    467468!-- 3D-array for storing the dissipation, needed for calculating the sgs
    468469!-- particle velocities
    469     IF ( use_sgs_for_particles  .OR.  wang_kernel )  THEN
     470    IF ( use_sgs_for_particles  .OR.  wang_kernel  .OR.  turbulence )  THEN
    470471       ALLOCATE ( diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    471472    ELSE
  • palm/trunk/SOURCE/init_cloud_physics.f90

    r1054 r1065  
    2020! Current revisions:
    2121! -----------------
     22! The Courant number of sedimentation can be controlled with c_sedimentation.
    2223!
    2324! Former revisions:
     
    8889!-- Calculate timestep according to precipitation
    8990    IF ( icloud_scheme == 0  .AND.  precipitation )  THEN
    90        dt_precipitation = MINVAL( dzu(nzb+2:nzt) ) / w_precipitation
     91       dt_precipitation = c_sedimentation * MINVAL( dzu(nzb+2:nzt) ) /        &
     92                          w_precipitation
    9193    ENDIF
    9294!
  • palm/trunk/SOURCE/microphysics.f90

    • Property svn:keywords set to Id
    r1054 r1065  
    44! Current revisions:
    55! -----------------
     6! Sedimentation process implemented according to Stevens and Seifert (2008).
     7! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens
     8! and Stevens, 2010).
    69!
    710! Former revisions:
     
    241244
    242245       INTEGER ::  i, j, k
    243        REAL    :: dr_min = 2.0E-6, dr_max = 1.0E-3
    244246
    245247       DO  k = nzb_2d(j,i)+1, nzt
    246           IF ( ( qr(k,j,i) > eps_sb ) .AND. ( nr(k,j,i) > eps_sb ) )  THEN
     248
     249          IF ( qr(k,j,i) <= eps_sb )  THEN
     250             qr(k,j,i) = 0.0
     251          ELSE
     252!
     253!--          Adjust number of raindrops to avoid nonlinear effects in
     254!--          sedimentation and evaporation of rain drops due to too small or
     255!--          too big weights of rain drops (Stevens and Seifert, 2008).
     256             IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) )  THEN
     257                nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin
     258             ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) )  THEN
     259                nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax
     260             ENDIF
     261             xr(k) = hyrho(k) * qr(k,j,i) / nr(k,j,i)
    247262!
    248263!--          Weight averaged diameter of rain drops:
     
    250265                       dpirho_l )**( 1.0 / 3.0 )
    251266!
    252 !--          Adjust number of raindrops to avoid nonlinear effects in
    253 !--          sedimentation and evaporation of rain drops due to too small or
    254 !--          big diameters of rain drops (Stevens and Seifert, 2008).
    255              IF ( dr(k) < dr_min )  THEN
    256                 nr(k,j,i) = qr(k,j,i) * hyrho(k) / dr_min**3 * dpirho_l
    257                 dr(k)     = dr_min
    258              ELSEIF ( dr(k) > dr_max )  THEN
    259                 nr(k,j,i) = qr(k,j,i) * hyrho(k) / dr_max**3 * dpirho_l
    260                 dr(k)     = dr_max
    261              ENDIF
    262 !
    263 !--          Mean weight of rain drops (Seifert and Beheng, 2006):
    264              xr(k) = MIN( MAX( hyrho(k) * qr(k,j,i) / nr(k,j,i), xrmin ), xrmax)
    265 !
    266267!--          Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
    267268!--          Stevens and Seifert, 2008):
    268              IF ( .NOT. mu_constant )  THEN
    269                 mu_r(k) = 10.0 * ( 1.0 + TANH( 1.2E3 * ( dr(k) - 1.4E-3 ) ) )
    270              ELSE
    271                 mu_r(k) = mu_constant_value
    272              ENDIF
     269             mu_r(k) = 10.0 * ( 1.0 + TANH( 1.2E3 * ( dr(k) - 1.4E-3 ) ) )
    273270!
    274271!--          Slope parameter of gamma distribution (Seifert, 2008):
     
    288285       USE control_parameters
    289286       USE statistics
     287       USE grid_variables
    290288   
    291289       IMPLICIT NONE
    292290
    293291       INTEGER ::  i, j, k
    294        REAL    ::  k_au, autocon, phi_au, tau_cloud, xc, nu_c
     292       REAL    ::  k_au, autocon, phi_au, tau_cloud, xc, nu_c, rc,   &
     293                   l_mix, re_lambda, alpha_cc, r_cc, sigma_cc, epsilon
    295294
    296295       k_au = k_cc / ( 20.0 * x0 )
     
    310309!--          Shape parameter of gamma distribution (Geoffroy et al., 2010):
    311310!--          (Use constant nu_c = 1.0 instead?)
    312              nu_c      = 1580.0 * hyrho(k) * ql(k,j,i) - 0.28
     311             nu_c      = 1.0 !MAX( 0.0, 1580.0 * hyrho(k) * ql(k,j,i) - 0.28 )
    313312!
    314313!--          Mean weight of cloud droplets:
    315314             xc        = hyrho(k) * ql(k,j,i) / nc
    316315!
     316!--          Parameterized turbulence effects on autoconversion (Seifert,
     317!--          Nuijens and Stevens, 2010)
     318             IF ( turbulence )  THEN
     319!
     320!--             Weight averaged radius of cloud droplets:
     321                rc = 0.5 * ( xc * dpirho_l )**( 1.0 / 3.0 )
     322
     323                alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0 + a_3 * nu_c )
     324                r_cc     = ( b_1 + b_2 * nu_c ) / ( 1.0 + b_3 * nu_c )
     325                sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0 + c_3 * nu_c )
     326!
     327!--             Mixing length (neglecting distance to ground and stratification)
     328                l_mix = ( dx * dy * dzu(k) )**( 1.0 / 3.0 )
     329!
     330!--             Limit dissipation rate according to Seifert, Nuijens and
     331!--             Stevens (2010)
     332                epsilon = MIN( 0.06, diss(k,j,i) )
     333!
     334!--             Compute Taylor-microscale Reynolds number:
     335                re_lambda = 6.0 / 11.0 * ( l_mix / c_const )**( 2.0 / 3.0 ) *  &
     336                            SQRT( 15.0 / kin_vis_air ) * epsilon**( 1.0 / 6.0 )
     337!
     338!--             The factor of 1.0E4 is needed to convert the dissipation rate
     339!--             from m2 s-3 to cm2 s-3.
     340                k_au = k_au * ( 1.0 +                                          &
     341                       epsilon * 1.0E4 * ( re_lambda * 1.0E-3 )**0.25 *        &
     342                       ( alpha_cc * EXP( -1.0 * ( ( rc - r_cc ) /              &
     343                       sigma_cc )**2 ) + beta_cc ) )
     344             ENDIF
     345!
    317346!--          Autoconversion rate (Seifert and Beheng, 2006):
    318              autocon   = k_au * ( nu_c + 2.0 ) * ( nu_c + 4.0 ) /             &
    319                          ( nu_c + 1.0 )**2 * ql(k,j,i)**2 * xc**2 *           &
    320                          ( 1.0 + phi_au / ( 1.0 - tau_cloud )**2 ) *          &
     347             autocon   = k_au * ( nu_c + 2.0 ) * ( nu_c + 4.0 ) /              &
     348                         ( nu_c + 1.0 )**2 * ql(k,j,i)**2 * xc**2 *            &
     349                         ( 1.0 + phi_au / ( 1.0 - tau_cloud )**2 ) *           &
    321350                         rho_surface
    322              autocon   = MIN( autocon, ql(k,j,i) / ( dt_3d *                  &
     351             autocon   = MIN( autocon, ql(k,j,i) / ( dt_3d *                   &
    323352                              weight_substep(intermediate_timestep_count) ) )
    324353!
     
    346375
    347376       INTEGER ::  i, j, k
    348        REAL    ::  accr, phi_ac, tau_cloud
     377       REAL    ::  accr, phi_ac, tau_cloud, k_cr
    349378
    350379       DO  k = nzb_2d(j,i)+1, nzt
     
    356385!--          Universal function for accretion process
    357386!--          (Seifert and Beheng, 2001):
    358              phi_ac    = tau_cloud / ( tau_cloud + 5.0E-5 )
    359              phi_ac    = ( phi_ac**2 )**2
     387             phi_ac = tau_cloud / ( tau_cloud + 5.0E-5 )
     388             phi_ac = ( phi_ac**2 )**2
     389!
     390!--          Parameterized turbulence effects on autoconversion (Seifert,
     391!--          Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to
     392!--          convert the dissipation (diss) from m2 s-3 to cm2 s-3.
     393             IF ( turbulence )  THEN
     394                k_cr = k_cr0 * ( 1.0 + 0.05 *                            &
     395                                 MIN( 600.0, diss(k,j,i) * 1.0E4 )**0.25 )
     396             ELSE
     397                k_cr = k_cr0                       
     398             ENDIF
    360399!
    361400!--          Accretion rate (Seifert and Beheng, 2006):
    362              accr      = k_cr * ql(k,j,i) * qr(k,j,i) * phi_ac *              &
    363                          SQRT( rho_surface * hyrho(k) )
    364              accr      = MIN( accr, ql(k,j,i) / ( dt_3d *                     &
    365                               weight_substep(intermediate_timestep_count) ) )
     401             accr = k_cr * ql(k,j,i) * qr(k,j,i) * phi_ac *              &
     402                    SQRT( rho_surface * hyrho(k) )
     403             accr = MIN( accr, ql(k,j,i) / ( dt_3d *                     &
     404                    weight_substep(intermediate_timestep_count) ) )
    366405!
    367406!--          Tendencies for q, qr, pt:
     
    390429
    391430       DO  k = nzb_2d(j,i)+1, nzt
    392           IF ( ( qr(k,j,i) > eps_sb )  .AND.  ( nr(k,j,i) > eps_sb ) )  THEN
     431          IF ( qr(k,j,i) > eps_sb )  THEN
    393432!
    394433!--          Selfcollection rate (Seifert and Beheng, 2006):
    395434!--          pirho_l**( 1.0 / 3.0 ) is necessary to convert [lambda_r] = m-1 to
    396435!--          kg**( 1.0 / 3.0 ).
    397              phi_sc = ( 1.0 + kappa_rr / lambda_r(k) *                        &
    398                       pirho_l**( 1.0 / 3.0 ) )**( -9 )
     436             phi_sc = 1.0 !( 1.0 + kappa_rr / lambda_r(k) *                        &
     437                      !pirho_l**( 1.0 / 3.0 ) )**( -9 )
    399438
    400439             selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * phi_sc *               &
     
    438477
    439478       DO  k = nzb_2d(j,i)+1, nzt
    440           IF ( ( qr(k,j,i) > eps_sb )  .AND.  ( nr(k,j,i) > eps_sb ) )  THEN
     479          IF ( qr(k,j,i) > eps_sb )  THEN
    441480!
    442481!--          Actual liquid water temperature:
     
    555594       IMPLICIT NONE
    556595
    557        INTEGER ::  i, j, k, n, n_substep
    558        REAL    ::  sed_nr_tend, sed_qr_tend
     596       INTEGER ::  i, j, k, k_run, n, n_substep
     597       REAL    ::  c_run, d_max, d_mean, d_min, dt_sedi, flux, mean, z_run
     598
     599       REAL, DIMENSION(nzb:nzt) :: c_nr, c_qr, d_nr, d_qr, nr_slope, qr_slope, &
     600                                   w_nr, w_qr
     601!
     602!--    Computation of sedimentation flux. Implementation according to Stevens
     603!--    and Seifert (2008).
    559604
    560605       IF ( intermediate_timestep_count == 1 )  prr(:,j,i) = 0.0
    561606
    562           sed_nr = 0.0
    563           sed_qr = 0.0
    564 
    565           DO  k = nzb_2d(j,i)+1, nzt
    566              IF ( ( qr(k,j,i) > eps_sb )  .AND.  ( nr(k,j,i) > eps_sb ) )  THEN
    567 !
    568 !--             Sedimentation of rain water content and rain drop concentration
    569 !--             according to Stevens and Seifert (2008):
    570                 sed_nr(k) = MIN( 20.0, MAX( 0.1, a_term - b_term * ( 1.0 +    &
    571                             c_term / lambda_r(k) )**( -1.0 *                  &
    572                             ( mu_r(k) + 1.0 ) ) ) ) * nr(k,j,i)
    573                 sed_qr(k) = MIN( 20.0, MAX( 0.1, a_term - b_term * ( 1.0 +    &
    574                             c_term / lambda_r(k) )**( -1.0 *                  &
    575                             ( mu_r(k) + 4.0 ) ) ) ) * qr(k,j,i) * hyrho(k)
    576 !
    577 !--             Computation of rain rate
    578                 prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) *              &
    579                              weight_substep(intermediate_timestep_count)
    580              ENDIF
    581           ENDDO
    582        
    583           DO  k = nzb_2d(j,i)+1, nzt
    584              sed_nr_tend = MAX( ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1),      &
    585                                 -nr(k,j,i) / ( dt_3d *                        &
    586                                 weight_substep(intermediate_timestep_count) ) )
    587              sed_qr_tend = MAX( ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) /     &
    588                                 hyrho(k),                                     &
    589                                 -qr(k,j,i) / ( dt_3d *                        &
    590                                 weight_substep(intermediate_timestep_count) ) )
    591 
    592              tend_nr(k,j,i) = tend_nr(k,j,i) + sed_nr_tend
    593              tend_qr(k,j,i) = tend_qr(k,j,i) + sed_qr_tend 
    594           ENDDO
     607       dt_sedi = dt_3d * weight_substep(intermediate_timestep_count)
     608
     609       w_nr = 0.0
     610       w_qr = 0.0
     611!
     612!--    Compute velocities
     613       DO  k = nzb_s_inner(j,i)+1, nzt
     614          IF ( qr(k,j,i) > eps_sb )  THEN
     615             w_nr(k) = MAX( 0.1, MIN( 20.0, a_term - b_term * ( 1.0 +          &
     616                       c_term / lambda_r(k) )**( -1.0 * ( mu_r(k) + 1.0 ) ) ) )
     617             w_qr(k) = MAX( 0.1, MIN( 20.0, a_term - b_term * ( 1.0 +          &
     618                       c_term / lambda_r(k) )**( -1.0 * ( mu_r(k) + 4.0 ) ) ) )
     619          ELSE
     620             w_nr(k) = 0.0
     621             w_qr(k) = 0.0
     622          ENDIF
     623       ENDDO
     624!
     625!--    Adjust boundary values
     626       w_nr(nzb_2d(j,i)) = w_nr(nzb_2d(j,i)+1)
     627       w_qr(nzb_2d(j,i)) = w_qr(nzb_2d(j,i)+1)
     628       w_nr(nzt) = w_nr(nzt-1)
     629       w_qr(nzt) = w_qr(nzt-1)
     630!
     631!--    Compute Courant number
     632       DO  k = nzb_s_inner(j,i)+1, nzt-1
     633          c_nr(k) = 0.25 * ( w_nr(k-1) + 2.0 * w_nr(k) + w_nr(k+1) ) * &
     634                    dt_sedi * ddzu(k)
     635          c_qr(k) = 0.25 * ( w_qr(k-1) + 2.0 * w_qr(k) + w_qr(k+1) ) * &
     636                    dt_sedi * ddzu(k)
     637       ENDDO
     638!
     639!--    Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):
     640       IF ( limiter_sedimentation )  THEN
     641
     642          qr(nzb_s_inner(j,i),j,i) = 0.0
     643          nr(nzb_s_inner(j,i),j,i) = 0.0
     644          qr(nzt,j,i) = 0.0
     645          nr(nzt,j,i) = 0.0
     646
     647          DO k = nzb_s_inner(j,i)+1, nzt-1
     648             d_mean = 0.5 * ( qr(k+1,j,i) + qr(k-1,j,i) )
     649             d_min  = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )
     650             d_max  = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)
     651
     652             qr_slope(k) = SIGN(1.0, d_mean) * MIN ( 2.0 * d_min, 2.0 * d_max, &
     653                                                     ABS( d_mean ) )
     654
     655             d_mean = 0.5 * ( nr(k+1,j,i) + nr(k-1,j,i) )
     656             d_min  = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )
     657             d_max  = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)
     658
     659             nr_slope(k) = SIGN(1.0, d_mean) * MIN ( 2.0 * d_min, 2.0 * d_max, &
     660                                                     ABS( d_mean ) )
     661          ENDDO
     662
     663       ELSE
     664          nr_slope = 0.0
     665          qr_slope = 0.0
     666       ENDIF
     667!
     668!--    Compute sedimentation flux
     669       DO  k = nzt-2, nzb_s_inner(j,i)+1, -1
     670!
     671!--       Sum up all rain drop number densities which contribute to the flux
     672!--       through k-1/2
     673          flux  = 0.0
     674          z_run = 0.0 ! height above z(k)
     675          k_run = k
     676          c_run = MIN( 1.0, c_nr(k) )
     677          DO WHILE ( c_run > 0.0  .AND.  k_run <= nzt-1 )
     678             flux  = flux + hyrho(k_run) *                                    &
     679                     ( nr(k_run,j,i) + nr_slope(k_run) * ( 1.0 - c_run ) *    &
     680                     0.5 ) * c_run * dzu(k_run)
     681             z_run = z_run + dzu(k_run)
     682             k_run = k_run + 1
     683             c_run = MIN( 1.0, c_nr(k_run) - z_run * ddzu(k_run) )
     684          ENDDO
     685!
     686!--       It is not allowed to sediment more rain drop number density than
     687!--       available
     688          flux = MIN( flux,                                                   &
     689                      hyrho(k) * dzu(k) * nr(k,j,i) + sed_nr(k+1) * dt_sedi )
     690
     691          sed_nr(k)      = flux / dt_sedi
     692          tend_nr(k,j,i) = tend_nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) *     &
     693                                            ddzu(k+1) / hyrho(k)
     694!
     695!--       Sum up all rain water content which contributes to the flux
     696!--       through k-1/2
     697          flux  = 0.0
     698          z_run = 0.0 ! height above z(k)
     699          k_run = k
     700          c_run = MIN( 1.0, c_qr(k) )
     701          DO WHILE ( c_run > 0.0  .AND.  k_run <= nzt-1 )
     702             flux  = flux + hyrho(k_run) *                                    &
     703                     ( qr(k_run,j,i) + qr_slope(k_run) * ( 1.0 - c_run ) *    &
     704                     0.5 ) * c_run * dzu(k_run)
     705             z_run = z_run + dzu(k_run)
     706             k_run = k_run + 1
     707             c_run = MIN( 1.0, c_qr(k_run) - z_run * ddzu(k_run) )
     708          ENDDO
     709!
     710!--       It is not allowed to sediment more rain water content than available
     711          flux = MIN( flux,                                                   &
     712                      hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) * dt_sedi )
     713
     714          sed_qr(k)      = flux / dt_sedi
     715          tend_qr(k,j,i) = tend_qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) *     &
     716                                            ddzu(k+1) / hyrho(k)
     717!
     718!--       Compute the rain rate
     719          prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) *                    &
     720                       weight_substep(intermediate_timestep_count)
     721       ENDDO
    595722!
    596723!--    Precipitation amount
     
    616743       IMPLICIT NONE
    617744 
    618        REAL    ::  gamm, xx,                                        &
    619                    ser, tmp, x_gamm, y_gamm
     745       REAL    ::  gamm, ser, tmp, x_gamm, xx, y_gamm
    620746       INTEGER ::  j
    621747 
  • palm/trunk/SOURCE/modules.f90

    r1054 r1065  
    2020! Current revisions:
    2121! -----------------
     22! + c_sedimentation, limiter_sedimentation, turbulence, a_1, a_2, a_3, b_1, b_2,
     23! + b_3, c_1, c_2, c_3, beta_cc
     24!
     25! bottom boundary condition of qr, nr changed from Dirichlet to Neumann
    2226!
    2327! Former revisions:
     
    468472
    469473    LOGICAL ::  curvature_solution_effects = .FALSE., &
    470                 ventilation_effect = .FALSE., &
    471                 mu_constant = .FALSE.
    472 
    473     REAL ::  a_vent = 0.78,      & ! coef. for ventilation effect
    474              a_term = 9.65,      & ! coef. for terminal velocity (m s-1)
    475              b_vent = 0.308,     & ! coef. for ventilation effect
    476              b_term = 9.8,       & ! coef. for terminal velocity (m s-1)
    477              bfactor,            &
    478              c_evap = 0.7,       & ! constant in evaporation
    479              c_term = 600.0,     & ! coef. for terminal velocity (m-1)
     474                limiter_sedimentation = .TRUE., &
     475                ventilation_effect = .FALSE.
     476               
     477
     478    REAL ::  a_1 = 8.69E-4,     & ! coef. in turb. parametrization (cm-2 s3)
     479             a_2 = -7.38E-5,    & ! coef. in turb. parametrization (cm-2 s3)
     480             a_3 = -1.40E-2,    & ! coef. in turb. parametrization
     481             a_term = 9.65,     & ! coef. for terminal velocity (m s-1)
     482             a_vent = 0.78,     & ! coef. for ventilation effect
     483             b_1 = 11.45E-6,    & ! coef. in turb. parametrization (m)
     484             b_2 = 9.68E-6,     & ! coef. in turb. parametrization (m)
     485             b_3 = 0.62,        & ! coef. in turb. parametrization
     486             b_term = 9.8,      & ! coef. for terminal velocity (m s-1)
     487             b_vent = 0.308,    & ! coef. for ventilation effect
     488             beta_cc = 3.09E-4, & ! coef. in turb. parametrization (cm-2 s3)
     489             bfactor,           &
     490             c_1 = 4.82E-6,     & ! coef. in turb. parametrization (m)
     491             c_2 = 4.8E-6,      & ! coef. in turb. parametrization (m)
     492             c_3 = 0.76,        & ! coef. in turb. parametrization
     493             c_const = 0.93,    & ! const. in Taylor-microscale Reynolds number
     494             c_evap = 0.7,      & ! constant in evaporation
     495             c_sedimentation = 2.0, & ! Courant number of sedimentation process
     496             c_term = 600.0,    & ! coef. for terminal velocity (m-1)
    480497             cof(6) = (/ 76.18009172947146,      & ! coefficients in the
    481498                         -86.50532032941677,     & ! numerical
     
    484501                         0.1208650973866179E-2,  &
    485502                         -0.5395239384953E-5 /), &
    486              cp = 1005.0,        & ! heat capacity of dry air (J kg-1 K-1)
     503             cp = 1005.0,       & ! heat capacity of dry air (J kg-1 K-1)
    487504             diff_coeff_l = 0.23E-4, & ! diffusivity of water vapor (m2 s-1)
    488505             effective_coll_efficiency, &
    489              eps_ros = 1.0E-4,   & ! accuracy of Rosenbrock method
    490              eps_sb = 1.0E-20,   & ! threshold in two-moments scheme
    491              k_cc = 4.44E09,     & ! const. rain-rain kernel (m3 kg-2 s-1)
    492              k_cr = 5.25,        & ! const. cloud-rain kernel (m3 kg-1 s-1)
    493              k_rr = 7.12,        & ! const. rain-rain kernel (m3 kg-1 s-1)
    494              k_br = 1000.,       & ! const. in breakup parametrization (m-1)
    495              kappa_rr = 60.7,    & ! const. in collision kernel (kg-1/3)
     506             eps_ros = 1.0E-4,  & ! accuracy of Rosenbrock method
     507             eps_sb = 1.0E-20,  & ! threshold in two-moments scheme
     508             k_cc = 9.44E09,    & ! const. rain-rain kernel (m3 kg-2 s-1)
     509             k_cr0 = 4.33,      & ! const. cloud-rain kernel (m3 kg-1 s-1)
     510             k_rr = 7.12,       & ! const. rain-rain kernel (m3 kg-1 s-1)
     511             k_br = 1000.,      & ! const. in breakup parametrization (m-1)
     512             kappa_rr = 60.7,   & ! const. in collision kernel (kg-1/3)
    496513             kin_vis_air = 1.4086E-5, & ! kin. viscosity of air (m2 s-1)
    497              l_v = 2.5E+06,      & ! latent heat of vaporization (J kg-1)
     514             l_v = 2.5E+06,     & ! latent heat of vaporization (J kg-1)
    498515             l_d_cp, l_d_r, l_d_rv, & ! l_v / cp, l_v / r_d, l_v / r_v
    499516             mass_of_solute = 1.0E-17, & ! soluted NaCl (kg)
    500517             molecular_weight_of_solute = 0.05844, & ! mol. m. NaCl (kg mol-1)
    501518             molecular_weight_of_water = 0.01801528, & ! mol. m. H2O (kg mol-1)
    502              mu_constant_value = 0.0, & ! shape param. of gamma distribution
    503              nc = 70.0E6,        & ! cloud droplet concentration
     519             nc = 70.0E6,       & ! cloud droplet concentration
    504520             prec_time_const = 0.001, & !coef. in Kessler scheme
    505              pirho_l, dpirho_l,  & ! pi * rho_l / 6.0; 6.0 / ( pi * rho_l )
    506              rho_l = 1.0E3,      & ! density of water (kg m-3)
    507              ql_crit = 0.0005,   & ! coef. in Kessler scheme
    508              r_d = 287.0,        & ! sp. gas const. dry air (J kg-1 K-1)
    509              r_v = 461.51,       & ! sp. gas const. water vapor (J kg-1 K-1)
    510              schmidt = 0.71,     & ! Schmidt number
    511              schmidt_p_1d3,      & ! schmidt**( 1.0 / 3.0 )
     521             pirho_l, dpirho_l, & ! pi * rho_l / 6.0; 6.0 / ( pi * rho_l )
     522             rho_l = 1.0E3,     & ! density of water (kg m-3)
     523             ql_crit = 0.0005,  & ! coef. in Kessler scheme
     524             r_d = 287.0,       & ! sp. gas const. dry air (J kg-1 K-1)
     525             r_v = 461.51,      & ! sp. gas const. water vapor (J kg-1 K-1)
     526             schmidt = 0.71,    & ! Schmidt number
     527             schmidt_p_1d3,     & ! schmidt**( 1.0 / 3.0 )
    512528             stp = 2.5066282746310005, & ! parameter in gamma function
    513529             thermal_conductivity_l = 2.43E-2, & ! therm. cond. air (J m-1 s-1 K-1)
    514              vanthoff = 2.0,     & ! van't Hoff factor for NaCl
    515              x0 = 2.6E-10,       & ! separating drop mass (kg)
    516              xrmin = 2.6E-10,    & ! minimum rain drop size (kg)
    517              xrmax = 5.0E-6,     & ! maximum rain drop site (kg)
     530             vanthoff = 2.0,    & ! van't Hoff factor for NaCl
     531             x0 = 2.6E-10,      & ! separating drop mass (kg)
     532             xrmin = 2.6E-10,   & ! minimum rain drop size (kg)
     533             xrmax = 5.0E-6,    & ! maximum rain drop site (kg)
    518534             dt_precipitation = 100.0, & ! timestep precipitation (s)
    519535             w_precipitation = 9.65      ! maximum terminal velocity (m s-1)
     
    595611                             scalar_advec = 'ws-scheme'
    596612    CHARACTER (LEN=20)   ::  bc_e_b = 'neumann', bc_lr = 'cyclic', &
    597                              bc_nr_b = 'dirichlet', bc_nr_t = 'neumann', &
     613                             bc_nr_b = 'neumann', bc_nr_t = 'neumann', &
    598614                             bc_ns = 'cyclic', bc_p_b = 'neumann', &
    599615                             bc_p_t = 'dirichlet', bc_pt_b = 'dirichlet', &
    600616                             bc_pt_t = 'initial_gradient', &
    601617                             bc_q_b = 'dirichlet', bc_q_t = 'neumann', &
    602                              bc_qr_b = 'dirichlet', bc_qr_t = 'neumann',&
     618                             bc_qr_b = 'neumann', bc_qr_t = 'neumann',&
    603619                             bc_s_b = 'dirichlet', bc_s_t = 'neumann', &
    604620                             bc_sa_t = 'neumann', &
     
    710726                do3d_compress = .FALSE., do_sum = .FALSE., &
    711727                dp_external = .FALSE., dp_smooth = .FALSE., &
    712                 drizzle = .TRUE., dt_fixed = .FALSE., &
     728                drizzle = .FALSE., dt_fixed = .FALSE., &
    713729                dt_3d_reached, dt_3d_reached_l, exchange_mg = .FALSE., &
    714730                first_call_lpm = .TRUE., &
     
    728744                scalar_rayleigh_damping = .TRUE., sloping_surface = .FALSE., &
    729745                stop_dt = .FALSE., synchronous_exchange = .FALSE., &
    730                 terminate_run = .FALSE., turbulent_inflow = .FALSE., &
     746                terminate_run = .FALSE., turbulence = .FALSE., &
     747                turbulent_inflow = .FALSE., &
    731748                use_prescribed_profile_data = .FALSE., use_reference = .FALSE.,&
    732749                use_surface_fluxes = .FALSE., use_top_fluxes = .FALSE., &
  • palm/trunk/SOURCE/parin.f90

    r1054 r1065  
    2020! Current revisions:
    2121! -----------------
     22! +nc, c_sedimentation, limiter_sedimentation, turbulence
     23! -mu_constant, mu_constant_value
    2224!
    2325! Former revisions:
     
    199201             call_psolver_at_all_substeps, canopy_mode, canyon_height, &
    200202             canyon_width_x, canyon_width_y, canyon_wall_left, &
    201              canyon_wall_south, cfl_factor, cloud_droplets, cloud_physics, &
    202              cloud_scheme, collective_wait, conserve_volume_flow, &
     203             canyon_wall_south, c_sedimentation, cfl_factor, cloud_droplets, &
     204             cloud_physics, cloud_scheme, collective_wait, &
     205             conserve_volume_flow, &
    203206             conserve_volume_flow_mode, coupling_start_time, cthf, &
    204207             curvature_solution_effects, cycle_mg, damp_level_1d, &
     
    212215             initializing_actions, km_constant, lad_surface, &
    213216             lad_vertical_gradient, lad_vertical_gradient_level, &
    214              leaf_surface_concentration, &
     217             leaf_surface_concentration, limiter_sedimentation, &
    215218             loop_optimization, masking_method, mg_cycles, &
    216219             mg_switch_to_pe0_level, mixing_length_1d, momentum_advec, &
    217              mu_constant, mu_constant_value, &
    218              netcdf_precision, neutral, ngsrb, nr_surface, &
     220             nc, netcdf_precision, neutral, ngsrb, nr_surface, &
    219221             nr_surface_initial_change, nr_vertical_gradient, &
    220222             nr_vertical_gradient_level, nsor, &
     
    241243             topography, topography_grid_convention, top_heatflux, &
    242244             top_momentumflux_u, top_momentumflux_v, top_salinityflux, &
    243              turbulent_inflow, ug_surface, ug_vertical_gradient, &
     245             turbulence, turbulent_inflow, ug_surface, ug_vertical_gradient, &
    244246             ug_vertical_gradient_level, use_surface_fluxes, &
    245247             use_top_fluxes, use_ug_for_galilei_tr, use_upstream_for_tke, &
  • palm/trunk/SOURCE/read_var_list.f90

    r1054 r1065  
    2020! Current revisions:
    2121! ------------------
     22! +nc, c_sedimentation, limiter_sedimentation, turbulence
     23! -mu_constant, mu_constant_value
    2224!
    2325! Former revisions:
     
    179181!-- Make version number check first
    180182    READ ( 13 )  version_on_file
    181     binary_version = '3.6'
     183    binary_version = '3.8'
    182184    IF ( TRIM( version_on_file ) /= TRIM( binary_version ) )  THEN
    183185       WRITE( message_string, * ) 'version mismatch concerning control ', &
     
    342344          CASE ( 'canyon_wall_south' )
    343345             READ ( 13 )  canyon_wall_south
     346          CASE ( 'c_sedimentation' )
     347             READ ( 13 )  c_sedimentation
    344348          CASE ( 'cfl_factor' )
    345349             READ ( 13 )  cfl_factor
     
    447451          CASE ( 'leaf_surface_concentration' )
    448452             READ ( 13 )  leaf_surface_concentration
     453          CASE ( 'limiter_sedimentation' )
     454             READ ( 13 )  limiter_sedimentation
    449455          CASE ( 'loop_optimization' )
    450456             READ ( 13 )  loop_optimization
     
    464470          CASE ( 'momentum_advec' )
    465471             READ ( 13 )  momentum_advec
    466           CASE ( 'mu_constant' )
    467              READ ( 13 )  mu_constant
    468           CASE ( 'mu_constant_value' )
    469              READ ( 13 )  mu_constant_value
     472          CASE ( 'nc' )
     473             READ ( 13 )  nc
    470474          CASE ( 'netcdf_precision' )
    471475             READ ( 13 )  netcdf_precision
     
    678682          CASE ( 'tsc' )
    679683             READ ( 13 )  tsc
     684          CASE ( 'turbulence' )
     685             READ ( 13 )  turbulence
    680686          CASE ( 'turbulent_inflow' )
    681687             READ ( 13 )  turbulent_inflow
  • palm/trunk/SOURCE/time_integration.f90

    r1054 r1065  
    2020! Current revisions:
    2121! -----------------
     22! exchange of diss (dissipation rate) in case of turbulence = .TRUE. added
    2223!
    2324! Former revisions:
     
    255256             CALL exchange_horiz( ql_vp, nbgp )
    256257          ENDIF
    257           IF ( wang_kernel )  CALL exchange_horiz( diss, nbgp )
     258          IF ( wang_kernel  .OR.  turbulence )  CALL exchange_horiz( diss, nbgp )
    258259
    259260          CALL cpu_log( log_point(26), 'exchange-horiz-progn', 'stop' )
  • palm/trunk/SOURCE/write_var_list.f90

    r1054 r1065  
    2020! Current revisions:
    2121! -----------------
     22! +nc, c_sedimentation, turbulence, limiter_sedimentation
     23! -mu_constant, mu_constant_value
    2224!
    2325! Former revisions:
     
    161163
    162164
    163     binary_version = '3.6'
     165    binary_version = '3.8'
    164166
    165167    WRITE ( 14 )  binary_version
     
    264266    WRITE ( 14 )  'canyon_wall_south             '
    265267    WRITE ( 14 )  canyon_wall_south
     268    WRITE ( 14 )  'c_sedimentation               '
     269    WRITE ( 14 )  c_sedimentation
    266270    WRITE ( 14 )  'cfl_factor                    '
    267271    WRITE ( 14 )  cfl_factor
     
    368372    WRITE ( 14 )  'leaf_surface_concentration    '
    369373    WRITE ( 14 )  leaf_surface_concentration
     374    WRITE ( 14 )  'limiter_sedimentation         '
     375    WRITE ( 14 )  limiter_sedimentation
    370376    WRITE ( 14 )  'loop_optimization             '
    371377    WRITE ( 14 )  loop_optimization
     
    384390    WRITE ( 14 )  'momentum_advec                '
    385391    WRITE ( 14 )  momentum_advec
    386     WRITE ( 14 )  'mu_constant                   '
    387     WRITE ( 14 )  mu_constant
    388     WRITE ( 14 )  'mu_constant_value             '
    389     WRITE ( 14 )  mu_constant_value
     392    WRITE ( 14 )  'nc                            '
     393    WRITE ( 14 )  nc
    390394    WRITE ( 14 )  'netcdf_precision              '
    391395    WRITE ( 14 )  netcdf_precision
     
    596600    WRITE ( 14 )  'tsc                           '
    597601    WRITE ( 14 )  tsc
     602    WRITE ( 14 )  'turbulence                    '
     603    WRITE ( 14 )  turbulence
    598604    WRITE ( 14 )  'turbulent_inflow              '
    599605    WRITE ( 14 )  turbulent_inflow
Note: See TracChangeset for help on using the changeset viewer.