Changeset 1353 for palm/trunk/SOURCE/microphysics.f90
- Timestamp:
- Apr 8, 2014 3:21:23 PM (10 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/microphysics.f90 (modified) (23 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/microphysics.f90
r1347 r1353 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! REAL constants provided with KIND-attribute 23 23 ! 24 24 ! Former revisions: … … 368 368 !-- t / pt : ratio of actual and potential temperature (t_d_pt) 369 369 !-- p_0(z) : vertical profile of the hydrostatic pressure (hyp) 370 t_surface = pt_surface * ( surface_pressure / 1000.0 )**0.286_wp370 t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp 371 371 DO k = nzb, nzt+1 372 hyp(k) = surface_pressure * 100.0 * &372 hyp(k) = surface_pressure * 100.0_wp * & 373 373 ( (t_surface - g/cp * zu(k)) / t_surface )**(1.0_wp/0.286_wp) 374 pt_d_t(k) = ( 100000.0 / hyp(k) )**0.286_wp375 t_d_pt(k) = 1.0 / pt_d_t(k)374 pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp 375 t_d_pt(k) = 1.0_wp / pt_d_t(k) 376 376 hyrho(k) = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) ) 377 377 ENDDO 378 378 ! 379 379 !-- Compute reference density 380 rho_surface = surface_pressure * 100.0 / ( r_d * t_surface )380 rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface ) 381 381 ENDIF 382 382 … … 445 445 446 446 IF ( qr(k,j,i) <= eps_sb ) THEN 447 qr(k,j,i) = 0.0 448 nr(k,j,i) = 0.0 447 qr(k,j,i) = 0.0_wp 448 nr(k,j,i) = 0.0_wp 449 449 ELSE 450 450 ! … … 506 506 REAL(wp) :: xc !: 507 507 508 k_au = k_cc / ( 20.0 * x0 )508 k_au = k_cc / ( 20.0_wp * x0 ) 509 509 510 510 DO k = nzb_s_inner(j,i)+1, nzt … … 513 513 ! 514 514 !-- Intern time scale of coagulation (Seifert and Beheng, 2006): 515 !-- (1.0 - qc(k,j,i) / ( qc(k,j,i) + qr_1d(k) ))516 tau_cloud = 1.0 - qc_1d(k) / ( qr_1d(k) + qc_1d(k) )515 !-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr_1d(k) )) 516 tau_cloud = 1.0_wp - qc_1d(k) / ( qr_1d(k) + qc_1d(k) ) 517 517 ! 518 518 !-- Universal function for autoconversion process 519 519 !-- (Seifert and Beheng, 2006): 520 phi_au = 600.0 * tau_cloud**0.68_wp * ( 1.0- tau_cloud**0.68_wp )**3520 phi_au = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**3 521 521 ! 522 522 !-- Shape parameter of gamma distribution (Geoffroy et al., 2010): 523 !-- (Use constant nu_c = 1.0 instead?)524 nu_c = 1.0 !MAX( 0.0_wp, 1580.0* hyrho(k) * qc(k,j,i) - 0.28_wp )523 !-- (Use constant nu_c = 1.0_wp instead?) 524 nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp ) 525 525 ! 526 526 !-- Mean weight of cloud droplets: … … 532 532 ! 533 533 !-- Weight averaged radius of cloud droplets: 534 rc = 0.5 * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )535 536 alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0 + a_3 * nu_c )537 r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0 + b_3 * nu_c )538 sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0 + c_3 * nu_c )534 rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp ) 535 536 alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c ) 537 r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c ) 538 sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c ) 539 539 ! 540 540 !-- Mixing length (neglecting distance to ground and stratification) … … 546 546 ! 547 547 !-- Compute Taylor-microscale Reynolds number: 548 re_lambda = 6.0 / 11.0* ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) * &549 SQRT( 15.0/ kin_vis_air ) * epsilon**( 1.0_wp / 6.0_wp )548 re_lambda = 6.0_wp / 11.0_wp * ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) * & 549 SQRT( 15.0_wp / kin_vis_air ) * epsilon**( 1.0_wp / 6.0_wp ) 550 550 ! 551 551 !-- The factor of 1.0E4 is needed to convert the dissipation rate 552 552 !-- from m2 s-3 to cm2 s-3. 553 k_au = k_au * ( 1.0_wp + &554 epsilon * 1.0E4 * ( re_lambda * 1.0E-3)**0.25_wp * &555 ( alpha_cc * EXP( -1.0 * ( ( rc - r_cc ) /&553 k_au = k_au * ( 1.0_wp + & 554 epsilon * 1.0E4_wp * ( re_lambda * 1.0E-3_wp )**0.25_wp * & 555 ( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) / & 556 556 sigma_cc )**2 ) + beta_cc ) ) 557 557 ENDIF 558 558 ! 559 559 !-- Autoconversion rate (Seifert and Beheng, 2006): 560 autocon = k_au * ( nu_c + 2.0 ) * ( nu_c + 4.0 ) /&561 ( nu_c + 1.0 )**2.0_wp * qc_1d(k)**2.0_wp * xc**2.0_wp * &562 ( 1.0 + phi_au / ( 1.0 - tau_cloud )**2.0_wp ) *&560 autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / & 561 ( nu_c + 1.0_wp )**2.0_wp * qc_1d(k)**2.0_wp * xc**2.0_wp * & 562 ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2.0_wp ) * & 563 563 rho_surface 564 564 autocon = MIN( autocon, qc_1d(k) / dt_micro ) … … 607 607 ! 608 608 !-- Intern time scale of coagulation (Seifert and Beheng, 2006): 609 tau_cloud = 1.0 - qc_1d(k) / ( qc_1d(k) + qr_1d(k) )609 tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) ) 610 610 ! 611 611 !-- Universal function for accretion process 612 612 !-- (Seifert and Beheng, 2001): 613 phi_ac = tau_cloud / ( tau_cloud + 5.0E-5 )613 phi_ac = tau_cloud / ( tau_cloud + 5.0E-5_wp ) 614 614 phi_ac = ( phi_ac**2.0_wp )**2.0_wp 615 615 ! … … 618 618 !-- convert the dissipation (diss) from m2 s-3 to cm2 s-3. 619 619 IF ( turbulence ) THEN 620 k_cr = k_cr0 * ( 1.0 + 0.05* &621 MIN( 600.0_wp, diss(k,j,i) * 1.0E4 )**0.25_wp )620 k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * & 621 MIN( 600.0_wp, diss(k,j,i) * 1.0E4_wp )**0.25_wp ) 622 622 ELSE 623 623 k_cr = k_cr0 … … 677 677 ! 678 678 !-- Collisional breakup rate (Seifert, 2008): 679 IF ( dr >= 0.3E-3 ) THEN680 phi_br = k_br * ( dr - 1.1E-3 )681 breakup = selfcoll * ( phi_br + 1.0 )679 IF ( dr >= 0.3E-3_wp ) THEN 680 phi_br = k_br * ( dr - 1.1E-3_wp ) 681 breakup = selfcoll * ( phi_br + 1.0_wp ) 682 682 ELSE 683 breakup = 0.0 683 breakup = 0.0_wp 684 684 ENDIF 685 685 … … 749 749 ! 750 750 !-- Saturation vapor pressure at t_l: 751 e_s = 610.78 * EXP( 17.269 * ( t_l - 273.16 ) / ( t_l - 35.86) )751 e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / ( t_l - 35.86_wp ) ) 752 752 ! 753 753 !-- Computation of saturation humidity: 754 q_s = 0.622 * e_s / ( hyp(k) - 0.378* e_s )755 alpha = 0.622 * l_d_r * l_d_cp / ( t_l * t_l )756 q_s = q_s * ( 1.0 + alpha * q_1d(k) ) / ( 1.0+ alpha * q_s )754 q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s ) 755 alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l ) 756 q_s = q_s * ( 1.0_wp + alpha * q_1d(k) ) / ( 1.0_wp + alpha * q_s ) 757 757 ! 758 758 !-- Supersaturation: 759 sat = MIN( 0.0_wp, ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0 )759 sat = MIN( 0.0_wp, ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp ) 760 760 ! 761 761 !-- Actual temperature: 762 762 temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) ) 763 763 764 g_evap = 1.0 / ( ( l_v / ( r_v * temp ) - 1.0) * l_v / &764 g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v / & 765 765 ( thermal_conductivity_l * temp ) + r_v * temp / & 766 766 ( diff_coeff_l * e_s ) ) … … 778 778 !-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005; 779 779 !-- Stevens and Seifert, 2008): 780 mu_r = 10.0 * ( 1.0 + TANH( 1.2E3 * ( dr - 1.4E-3) ) )780 mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) ) 781 781 ! 782 782 !-- Slope parameter of gamma distribution (Seifert, 2008): 783 lambda_r = ( ( mu_r + 3.0 ) * ( mu_r + 2.0) * &784 ( mu_r + 1.0 ) )**( 1.0_wp / 3.0_wp ) / dr785 786 mu_r_2 = mu_r + 2.0 787 mu_r_5d2 = mu_r + 2.5 783 lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * & 784 ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr 785 786 mu_r_2 = mu_r + 2.0_wp 787 mu_r_5d2 = mu_r + 2.5_wp 788 788 f_vent = a_vent * gamm( mu_r_2 ) * & 789 789 lambda_r**( -mu_r_2 ) + & … … 791 791 SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) * & 792 792 lambda_r**( -mu_r_5d2 ) * & 793 ( 1.0 - 0.5 * ( b_term / a_term ) *&793 ( 1.0_wp - 0.5_wp * ( b_term / a_term ) * & 794 794 ( lambda_r / & 795 795 ( c_term + lambda_r ) )**mu_r_5d2 - & 796 0.125 * ( b_term / a_term )**2.0_wp *&796 0.125_wp * ( b_term / a_term )**2.0_wp * & 797 797 ( lambda_r / & 798 ( 2.0 * c_term + lambda_r ) )**mu_r_5d2 -&799 0.0625 * ( b_term / a_term )**3.0_wp *&798 ( 2.0_wp * c_term + lambda_r ) )**mu_r_5d2 - & 799 0.0625_wp * ( b_term / a_term )**3.0_wp * & 800 800 ( lambda_r / & 801 ( 3.0 * c_term + lambda_r ) )**mu_r_5d2 -&802 0.0390625 * ( b_term / a_term )**4.0_wp *&801 ( 3.0_wp * c_term + lambda_r ) )**mu_r_5d2 - & 802 0.0390625_wp * ( b_term / a_term )**4.0_wp * & 803 803 ( lambda_r / & 804 ( 4.0 * c_term + lambda_r ) )**mu_r_5d2 )805 nr_0 = nr_1d(k) * lambda_r**( mu_r + 1.0 ) /&806 gamm( mu_r + 1.0 )804 ( 4.0_wp * c_term + lambda_r ) )**mu_r_5d2 ) 805 nr_0 = nr_1d(k) * lambda_r**( mu_r + 1.0_wp ) / & 806 gamm( mu_r + 1.0_wp ) 807 807 ELSE 808 f_vent = 1.0 808 f_vent = 1.0_wp 809 809 nr_0 = nr_1d(k) * dr 810 810 ENDIF 811 811 ! 812 812 !-- Evaporation rate of rain water content (Seifert and Beheng, 2006): 813 evap = 2.0 * pi * nr_0 * g_evap * f_vent * sat / &813 evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / & 814 814 hyrho(k) 815 815 … … 859 859 ! 860 860 !-- Sedimentation of cloud droplets (Heus et al., 2010): 861 sed_qc_const = k_st * ( 3.0 / ( 4.0* pi * rho_l ))**( 2.0_wp / 3.0_wp ) * &862 EXP( 5.0 * LOG( sigma_gc )**2 )863 864 sed_qc(nzt+1) = 0.0 861 sed_qc_const = k_st * ( 3.0_wp / ( 4.0_wp * pi * rho_l ))**( 2.0_wp / 3.0_wp ) * & 862 EXP( 5.0_wp * LOG( sigma_gc )**2 ) 863 864 sed_qc(nzt+1) = 0.0_wp 865 865 866 866 DO k = nzt, nzb_s_inner(j,i)+1, -1 … … 869 869 ( qc_1d(k) * hyrho(k) )**( 5.0_wp / 3.0_wp ) 870 870 ELSE 871 sed_qc(k) = 0.0 871 sed_qc(k) = 0.0_wp 872 872 ENDIF 873 873 … … 943 943 !-- Computation of sedimentation flux. Implementation according to Stevens 944 944 !-- and Seifert (2008). 945 IF ( intermediate_timestep_count == 1 ) prr(:,j,i) = 0.0 945 IF ( intermediate_timestep_count == 1 ) prr(:,j,i) = 0.0_wp 946 946 ! 947 947 !-- Compute velocities … … 954 954 !-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005; 955 955 !-- Stevens and Seifert, 2008): 956 mu_r = 10.0 * ( 1.0 + TANH( 1.2E3 * ( dr - 1.4E-3) ) )956 mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) ) 957 957 ! 958 958 !-- Slope parameter of gamma distribution (Seifert, 2008): 959 lambda_r = ( ( mu_r + 3.0 ) * ( mu_r + 2.0) * &960 ( mu_r + 1.0 ) )**( 1.0_wp / 3.0_wp ) / dr961 962 w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, a_term - b_term * ( 1.0 + &963 c_term / lambda_r )**( -1.0 * ( mu_r + 1.0) ) ) )964 w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, a_term - b_term * ( 1.0 + &965 c_term / lambda_r )**( -1.0 * ( mu_r + 4.0) ) ) )959 lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * & 960 ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr 961 962 w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, a_term - b_term * ( 1.0_wp + & 963 c_term / lambda_r )**( -1.0_wp * ( mu_r + 1.0_wp ) ) ) ) 964 w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, a_term - b_term * ( 1.0_wp + & 965 c_term / lambda_r )**( -1.0_wp * ( mu_r + 4.0_wp ) ) ) ) 966 966 ELSE 967 w_nr(k) = 0.0 968 w_qr(k) = 0.0 967 w_nr(k) = 0.0_wp 968 w_qr(k) = 0.0_wp 969 969 ENDIF 970 970 ENDDO … … 973 973 w_nr(nzb_s_inner(j,i)) = w_nr(nzb_s_inner(j,i)+1) 974 974 w_qr(nzb_s_inner(j,i)) = w_qr(nzb_s_inner(j,i)+1) 975 w_nr(nzt+1) = 0.0 976 w_qr(nzt+1) = 0.0 975 w_nr(nzt+1) = 0.0_wp 976 w_qr(nzt+1) = 0.0_wp 977 977 ! 978 978 !-- Compute Courant number 979 979 DO k = nzb_s_inner(j,i)+1, nzt 980 c_nr(k) = 0.25 * ( w_nr(k-1) + 2.0* w_nr(k) + w_nr(k+1) ) * &980 c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) * & 981 981 dt_micro * ddzu(k) 982 c_qr(k) = 0.25 * ( w_qr(k-1) + 2.0* w_qr(k) + w_qr(k+1) ) * &982 c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) * & 983 983 dt_micro * ddzu(k) 984 984 ENDDO … … 988 988 989 989 DO k = nzb_s_inner(j,i)+1, nzt 990 d_mean = 0.5 * ( qr_1d(k+1) + qr_1d(k-1) )990 d_mean = 0.5_wp * ( qr_1d(k+1) + qr_1d(k-1) ) 991 991 d_min = qr_1d(k) - MIN( qr_1d(k+1), qr_1d(k), qr_1d(k-1) ) 992 992 d_max = MAX( qr_1d(k+1), qr_1d(k), qr_1d(k-1) ) - qr_1d(k) 993 993 994 qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0 * d_min, 2.0* d_max, &994 qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, 2.0_wp * d_max, & 995 995 ABS( d_mean ) ) 996 996 997 d_mean = 0.5 * ( nr_1d(k+1) + nr_1d(k-1) )997 d_mean = 0.5_wp * ( nr_1d(k+1) + nr_1d(k-1) ) 998 998 d_min = nr_1d(k) - MIN( nr_1d(k+1), nr_1d(k), nr_1d(k-1) ) 999 999 d_max = MAX( nr_1d(k+1), nr_1d(k), nr_1d(k-1) ) - nr_1d(k) 1000 1000 1001 nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0 * d_min, 2.0* d_max, &1001 nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, 2.0_wp * d_max, & 1002 1002 ABS( d_mean ) ) 1003 1003 ENDDO … … 1005 1005 ELSE 1006 1006 1007 nr_slope = 0.0 1008 qr_slope = 0.0 1007 nr_slope = 0.0_wp 1008 qr_slope = 0.0_wp 1009 1009 1010 1010 ENDIF 1011 1011 1012 sed_nr(nzt+1) = 0.0 1013 sed_qr(nzt+1) = 0.0 1012 sed_nr(nzt+1) = 0.0_wp 1013 sed_qr(nzt+1) = 0.0_wp 1014 1014 ! 1015 1015 !-- Compute sedimentation flux … … 1018 1018 !-- Sum up all rain drop number densities which contribute to the flux 1019 1019 !-- through k-1/2 1020 flux = 0.0 1021 z_run = 0.0 ! height above z(k)1020 flux = 0.0_wp 1021 z_run = 0.0_wp ! height above z(k) 1022 1022 k_run = k 1023 1023 c_run = MIN( 1.0_wp, c_nr(k) ) 1024 DO WHILE ( c_run > 0.0 .AND. k_run <= nzt )1024 DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt ) 1025 1025 flux = flux + hyrho(k_run) * & 1026 ( nr_1d(k_run) + nr_slope(k_run) * ( 1.0 - c_run ) * &1027 0.5 ) * c_run * dzu(k_run)1026 ( nr_1d(k_run) + nr_slope(k_run) * ( 1.0_wp - c_run ) * & 1027 0.5_wp ) * c_run * dzu(k_run) 1028 1028 z_run = z_run + dzu(k_run) 1029 1029 k_run = k_run + 1 … … 1042 1042 !-- Sum up all rain water content which contributes to the flux 1043 1043 !-- through k-1/2 1044 flux = 0.0 1045 z_run = 0.0 ! height above z(k)1044 flux = 0.0_wp 1045 z_run = 0.0_wp ! height above z(k) 1046 1046 k_run = k 1047 1047 c_run = MIN( 1.0_wp, c_qr(k) ) 1048 1048 1049 DO WHILE ( c_run > 0.0 .AND. k_run <= nzt-1 )1049 DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt-1 ) 1050 1050 1051 1051 flux = flux + hyrho(k_run) * & 1052 ( qr_1d(k_run) + qr_slope(k_run) * ( 1.0 - c_run ) * &1053 0.5 ) * c_run * dzu(k_run)1052 ( qr_1d(k_run) + qr_slope(k_run) * ( 1.0_wp - c_run ) * & 1053 0.5_wp ) * c_run * dzu(k_run) 1054 1054 z_run = z_run + dzu(k_run) 1055 1055 k_run = k_run + 1 … … 1114 1114 x_gamm = xx 1115 1115 y_gamm = x_gamm 1116 tmp = x_gamm + 5.5 1117 tmp = ( x_gamm + 0.5 ) * LOG( tmp ) - tmp1116 tmp = x_gamm + 5.5_wp 1117 tmp = ( x_gamm + 0.5_wp ) * LOG( tmp ) - tmp 1118 1118 ser = 1.000000000190015_wp 1119 1119 1120 1120 DO j = 1, 6 1121 y_gamm = y_gamm + 1.0 1121 y_gamm = y_gamm + 1.0_wp 1122 1122 ser = ser + cof( j ) / y_gamm 1123 1123 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
