Changeset 1334 for palm/trunk/SOURCE


Ignore:
Timestamp:
Mar 25, 2014 12:21:40 PM (11 years ago)
Author:
heinze
Message:

Bugfix: REAL constants provided with KIND-attribute

Location:
palm/trunk/SOURCE
Files:
2 edited

Legend:

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

    r1323 r1334  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Bugfix: REAL constants provided with KIND-attribute
    2323!
    2424! Former revisions:
     
    118118!-- t / pt : ratio of actual and potential temperature (t_d_pt)
    119119!-- p_0(z) : vertical profile of the hydrostatic pressure (hyp)
    120     t_surface = pt_surface * ( surface_pressure / 1000.0 )**0.286
     120    t_surface = pt_surface * ( surface_pressure / 1000.0 )**0.286_wp
    121121    DO  k = nzb, nzt+1
    122122!
     
    129129       hyp(k)    = surface_pressure * 100.0 * &
    130130                   ( (t_surface - g/cp * zu(k)) / t_surface )**(1.0_wp/0.286_wp)
    131        pt_d_t(k) = ( 100000.0 / hyp(k) )**0.286
     131       pt_d_t(k) = ( 100000.0 / hyp(k) )**0.286_wp
    132132       t_d_pt(k) = 1.0 / pt_d_t(k)
    133133       hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )       
  • palm/trunk/SOURCE/microphysics.f90

    r1323 r1334  
    2020! Current revisions:
    2121! ------------------
    22 !
     22! Bugfix: REAL constants provided with KIND-attribute
    2323!
    2424! Former revisions:
     
    361361!--       t / pt : ratio of actual and potential temperature (t_d_pt)
    362362!--       p_0(z) : vertical profile of the hydrostatic pressure (hyp)
    363           t_surface = pt_surface * ( surface_pressure / 1000.0 )**0.286
     363          t_surface = pt_surface * ( surface_pressure / 1000.0 )**0.286_wp
    364364          DO  k = nzb, nzt+1
    365365             hyp(k)    = surface_pressure * 100.0 * &
    366                          ( (t_surface - g/cp * zu(k)) / t_surface )**(1.0/0.286)
    367              pt_d_t(k) = ( 100000.0 / hyp(k) )**0.286
     366                         ( (t_surface - g/cp * zu(k)) / t_surface )**(1.0_wp/0.286_wp)
     367             pt_d_t(k) = ( 100000.0 / hyp(k) )**0.286_wp
    368368             t_d_pt(k) = 1.0 / pt_d_t(k)
    369369             hyrho(k)  = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) )       
     
    511511!--          Universal function for autoconversion process
    512512!--          (Seifert and Beheng, 2006):
    513              phi_au    = 600.0 * tau_cloud**0.68 * ( 1.0 - tau_cloud**0.68 )**3
     513             phi_au    = 600.0 * tau_cloud**0.68_wp * ( 1.0 - tau_cloud**0.68_wp )**3
    514514!
    515515!--          Shape parameter of gamma distribution (Geoffroy et al., 2010):
    516516!--          (Use constant nu_c = 1.0 instead?)
    517              nu_c      = 1.0 !MAX( 0.0, 1580.0 * hyrho(k) * qc(k,j,i) - 0.28 )
     517             nu_c      = 1.0 !MAX( 0.0_wp, 1580.0 * hyrho(k) * qc(k,j,i) - 0.28_wp )
    518518!
    519519!--          Mean weight of cloud droplets:
     
    525525!
    526526!--             Weight averaged radius of cloud droplets:
    527                 rc = 0.5 * ( xc * dpirho_l )**( 1.0 / 3.0_wp )
     527                rc = 0.5 * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )
    528528
    529529                alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0 + a_3 * nu_c )
     
    532532!
    533533!--             Mixing length (neglecting distance to ground and stratification)
    534                 l_mix = ( dx * dy * dzu(k) )**( 1.0 / 3.0_wp )
     534                l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )
    535535!
    536536!--             Limit dissipation rate according to Seifert, Nuijens and
    537537!--             Stevens (2010)
    538                 epsilon = MIN( 0.06, diss(k,j,i) )
     538                epsilon = MIN( 0.06_wp, diss(k,j,i) )
    539539!
    540540!--             Compute Taylor-microscale Reynolds number:
    541                 re_lambda = 6.0 / 11.0 * ( l_mix / c_const )**( 2.0 / 3.0_wp ) *  &
    542                             SQRT( 15.0 / kin_vis_air ) * epsilon**( 1.0 / 6.0_wp )
     541                re_lambda = 6.0 / 11.0 * ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) *  &
     542                            SQRT( 15.0 / kin_vis_air ) * epsilon**( 1.0_wp / 6.0_wp )
    543543!
    544544!--             The factor of 1.0E4 is needed to convert the dissipation rate
    545545!--             from m2 s-3 to cm2 s-3.
    546                 k_au = k_au * ( 1.0 +                                          &
    547                        epsilon * 1.0E4 * ( re_lambda * 1.0E-3 )**0.25 *        &
     546                k_au = k_au * ( 1.0_wp +                                       &
     547                       epsilon * 1.0E4 * ( re_lambda * 1.0E-3 )**0.25_wp *     &
    548548                       ( alpha_cc * EXP( -1.0 * ( ( rc - r_cc ) /              &
    549549                       sigma_cc )**2 ) + beta_cc ) )
     
    552552!--          Autoconversion rate (Seifert and Beheng, 2006):
    553553             autocon = k_au * ( nu_c + 2.0 ) * ( nu_c + 4.0 ) /    &
    554                        ( nu_c + 1.0 )**2 * qc_1d(k)**2 * xc**2 *   &
    555                        ( 1.0 + phi_au / ( 1.0 - tau_cloud )**2 ) * &
     554                       ( nu_c + 1.0 )**2.0_wp * qc_1d(k)**2.0_wp * xc**2.0_wp *   &
     555                       ( 1.0 + phi_au / ( 1.0 - tau_cloud )**2.0_wp ) * &
    556556                       rho_surface
    557557             autocon = MIN( autocon, qc_1d(k) / dt_micro )
     
    605605!--          (Seifert and Beheng, 2001):
    606606             phi_ac = tau_cloud / ( tau_cloud + 5.0E-5 )
    607              phi_ac = ( phi_ac**2 )**2
     607             phi_ac = ( phi_ac**2.0_wp )**2.0_wp
    608608!
    609609!--          Parameterized turbulence effects on autoconversion (Seifert,
     
    612612             IF ( turbulence )  THEN
    613613                k_cr = k_cr0 * ( 1.0 + 0.05 *                             &
    614                                  MIN( 600.0, diss(k,j,i) * 1.0E4 )**0.25 )
     614                                 MIN( 600.0_wp, diss(k,j,i) * 1.0E4 )**0.25_wp )
    615615             ELSE
    616616                k_cr = k_cr0                       
     
    667667!
    668668!--          Weight averaged diameter of rain drops:
    669              dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0 / 3.0_wp )
     669             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
    670670!
    671671!--          Collisional breakup rate (Seifert, 2008):
     
    750750!
    751751!--          Supersaturation:
    752              sat = MIN( 0.0, ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0 )
     752             sat = MIN( 0.0_wp, ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0 )
    753753!
    754754!--          Actual temperature:
     
    763763!
    764764!--          Weight averaged diameter of rain drops:
    765              dr = ( xr * dpirho_l )**( 1.0 / 3.0_wp )
     765             dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp )
    766766!
    767767!--          Compute ventilation factor and intercept parameter
     
    775775!--             Slope parameter of gamma distribution (Seifert, 2008):
    776776                lambda_r = ( ( mu_r + 3.0 ) * ( mu_r + 2.0 ) *       &
    777                              ( mu_r + 1.0 ) )**( 1.0 / 3.0_wp ) / dr
     777                             ( mu_r + 1.0 ) )**( 1.0_wp / 3.0_wp ) / dr
    778778
    779779                mu_r_2   = mu_r + 2.0
     
    787787                         ( lambda_r /                                         &
    788788                         (       c_term + lambda_r ) )**mu_r_5d2 -            &
    789                                  0.125 * ( b_term / a_term )**2 *             &
     789                                 0.125 * ( b_term / a_term )**2.0_wp *        &
    790790                         ( lambda_r /                                         &
    791791                         ( 2.0 * c_term + lambda_r ) )**mu_r_5d2 -            &
    792                                  0.0625 * ( b_term / a_term )**3 *            &
     792                                 0.0625 * ( b_term / a_term )**3.0_wp *       &
    793793                         ( lambda_r /                                         &
    794794                         ( 3.0 * c_term + lambda_r ) )**mu_r_5d2 -            &
    795                                  0.0390625 * ( b_term / a_term )**4 *         &
     795                                 0.0390625 * ( b_term / a_term )**4.0_wp *    &
    796796                         ( lambda_r /                                         &
    797797                         ( 4.0 * c_term + lambda_r ) )**mu_r_5d2 )
     
    852852!
    853853!--    Sedimentation of cloud droplets (Heus et al., 2010):
    854        sed_qc_const = k_st * ( 3.0 / ( 4.0 * pi * rho_l ))**( 2.0 / 3.0_wp ) *   &
     854       sed_qc_const = k_st * ( 3.0 / ( 4.0 * pi * rho_l ))**( 2.0_wp / 3.0_wp ) *   &
    855855                     EXP( 5.0 * LOG( sigma_gc )**2 )
    856856
     
    859859       DO  k = nzt, nzb_s_inner(j,i)+1, -1
    860860          IF ( qc_1d(k) > eps_sb )  THEN
    861              sed_qc(k) = sed_qc_const * nc_1d(k)**( -2.0 / 3.0_wp ) * &
    862                         ( qc_1d(k) * hyrho(k) )**( 5.0 / 3.0_wp )
     861             sed_qc(k) = sed_qc_const * nc_1d(k)**( -2.0_wp / 3.0_wp ) * &
     862                        ( qc_1d(k) * hyrho(k) )**( 5.0_wp / 3.0_wp )
    863863          ELSE
    864864             sed_qc(k) = 0.0
     
    943943!
    944944!--          Weight averaged diameter of rain drops:
    945              dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0 / 3.0_wp )
     945             dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp )
    946946!
    947947!--          Shape parameter of gamma distribution (Milbrandt and Yau, 2005;
     
    951951!--          Slope parameter of gamma distribution (Seifert, 2008):
    952952             lambda_r = ( ( mu_r + 3.0 ) * ( mu_r + 2.0 ) *          &
    953                         ( mu_r + 1.0 ) )**( 1.0 / 3.0_wp ) / dr
    954 
    955              w_nr(k) = MAX( 0.1, MIN( 20.0, a_term - b_term * ( 1.0 +          &
     953                        ( mu_r + 1.0 ) )**( 1.0_wp / 3.0_wp ) / dr
     954
     955             w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, a_term - b_term * ( 1.0 +    &
    956956                       c_term / lambda_r )**( -1.0 * ( mu_r + 1.0 ) ) ) )
    957              w_qr(k) = MAX( 0.1, MIN( 20.0, a_term - b_term * ( 1.0 +          &
     957             w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, a_term - b_term * ( 1.0 +    &
    958958                       c_term / lambda_r )**( -1.0 * ( mu_r + 4.0 ) ) ) )
    959959          ELSE
     
    11091109       tmp = x_gamm + 5.5
    11101110       tmp = ( x_gamm + 0.5 ) * LOG( tmp ) - tmp
    1111        ser = 1.000000000190015 
     1111       ser = 1.000000000190015_wp
    11121112
    11131113       DO  j = 1, 6
Note: See TracChangeset for help on using the changeset viewer.