Changeset 2155 for palm/trunk/SOURCE/microphysics_mod.f90
- Timestamp:
- Feb 21, 2017 9:57:40 AM (7 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/microphysics_mod.f90 (modified) (85 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/microphysics_mod.f90
r2101 r2155 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 23 ! 22 ! Bugfix in the calculation of microphysical quantities on ghost points. 23 ! 24 24 ! Former revisions: 25 25 ! ----------------- … … 28 28 ! 2031 2016-10-21 15:11:58Z knoop 29 29 ! renamed variable rho to rho_ocean 30 ! 30 ! 31 31 ! 2000 2016-08-20 18:09:15Z knoop 32 32 ! Forced header and separation lines into 80 columns 33 ! 33 ! 34 34 ! 1850 2016-04-08 13:29:27Z maronga 35 35 ! Module renamed … … 52 52 ! Added new routine calc_precipitation_amount. The routine now allows to account 53 53 ! for precipitation due to sedimenation of cloud (fog) droplets 54 ! 54 ! 55 55 ! 1682 2015-10-07 23:56:08Z knoop 56 ! Code annotations made doxygen readable 56 ! Code annotations made doxygen readable 57 57 ! 58 58 ! 1646 2015-09-02 16:00:10Z hoffmann … … 63 63 ! Vectorized version of adjust_cloud added. 64 64 ! Little reformatting of the code. 65 ! 65 ! 66 66 ! 1353 2014-04-08 15:21:23Z heinze 67 67 ! REAL constants provided with KIND-attribute 68 ! 68 ! 69 69 ! 1346 2014-03-27 13:18:20Z heinze 70 ! Bugfix: REAL constants provided with KIND-attribute especially in call of 70 ! Bugfix: REAL constants provided with KIND-attribute especially in call of 71 71 ! intrinsic function like MAX, MIN, SIGN 72 ! 72 ! 73 73 ! 1334 2014-03-25 12:21:40Z heinze 74 74 ! Bugfix: REAL constants provided with KIND-attribute … … 106 106 ! 1053 2012-11-13 17:11:03Z hoffmann 107 107 ! initial revision 108 ! 108 ! 109 109 ! Description: 110 110 ! ------------ … … 140 140 REAL(wp) :: c_term = 600.0_wp !< coef. for terminal velocity (m-1) 141 141 REAL(wp) :: diff_coeff_l = 0.23E-4_wp !< diffusivity of water vapor (m2 s-1) 142 REAL(wp) :: eps_sb = 1.0E- 20_wp !< threshold in two-moments scheme142 REAL(wp) :: eps_sb = 1.0E-10_wp !< threshold in two-moments scheme 143 143 REAL(wp) :: k_cc = 9.44E09_wp !< const. cloud-cloud kernel (m3 kg-2 s-1) 144 144 REAL(wp) :: k_cr0 = 4.33_wp !< const. cloud-rain kernel (m3 kg-1 s-1) … … 221 221 MODULE PROCEDURE sedimentation_cloud_ij 222 222 END INTERFACE sedimentation_cloud 223 223 224 224 INTERFACE sedimentation_rain 225 225 MODULE PROCEDURE sedimentation_rain … … 337 337 pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp 338 338 t_d_pt(k) = 1.0_wp / pt_d_t(k) 339 hyrho(k) = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) ) 339 hyrho(k) = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) ) 340 340 ENDDO 341 341 … … 346 346 347 347 ! 348 !-- Compute length of time step 348 !-- Compute length of time step 349 349 IF ( call_microphysics_at_all_substeps ) THEN 350 350 dt_micro = dt_3d * weight_pres(intermediate_timestep_count) … … 383 383 ! Description: 384 384 ! ------------ 385 !> Adjust number of raindrops to avoid nonlinear effects in sedimentation and 386 !> evaporation of rain drops due to too small or too big weights 385 !> Adjust number of raindrops to avoid nonlinear effects in sedimentation and 386 !> evaporation of rain drops due to too small or too big weights 387 387 !> of rain drops (Stevens and Seifert, 2008). 388 388 !------------------------------------------------------------------------------! … … 399 399 400 400 USE indices, & 401 ONLY: nxl , nxr, nys, nyn, nzb_s_inner, nzt401 ONLY: nxlg, nxrg, nysg, nyng, nzb_s_inner, nzt 402 402 403 403 USE kinds … … 411 411 CALL cpu_log( log_point_s(54), 'adjust_cloud', 'start' ) 412 412 413 DO i = nxl , nxr414 DO j = nys , nyn413 DO i = nxlg, nxrg 414 DO j = nysg, nyng 415 415 DO k = nzb_s_inner(j,i)+1, nzt 416 416 IF ( qr(k,j,i) <= eps_sb ) THEN … … 456 456 457 457 USE indices, & 458 ONLY: nxl , nxr, nys, nyn, nzb_s_inner, nzt458 ONLY: nxlg, nxrg, nysg, nyng, nzb_s_inner, nzt 459 459 460 460 USE kinds … … 466 466 INTEGER(iwp) :: k !< 467 467 468 REAL(wp) :: alpha_cc !< 468 REAL(wp) :: alpha_cc !< 469 469 REAL(wp) :: autocon !< 470 470 REAL(wp) :: dissipation !< … … 482 482 CALL cpu_log( log_point_s(55), 'autoconversion', 'start' ) 483 483 484 DO i = nxl , nxr485 DO j = nys , nyn484 DO i = nxlg, nxrg 485 DO j = nysg, nyng 486 486 DO k = nzb_s_inner(j,i)+1, nzt 487 487 … … 494 494 tau_cloud = 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) ) 495 495 ! 496 !-- Universal function for autoconversion process 496 !-- Universal function for autoconversion process 497 497 !-- (Seifert and Beheng, 2006): 498 498 phi_au = 600.0_wp * tau_cloud**0.68_wp * & … … 506 506 xc = hyrho(k) * qc(k,j,i) / nc_const 507 507 ! 508 !-- Parameterized turbulence effects on autoconversion (Seifert, 508 !-- Parameterized turbulence effects on autoconversion (Seifert, 509 509 !-- Nuijens and Stevens, 2010) 510 510 IF ( collision_turbulence ) THEN … … 517 517 sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c ) 518 518 ! 519 !-- Mixing length (neglecting distance to ground and 519 !-- Mixing length (neglecting distance to ground and 520 520 !-- stratification) 521 521 l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp ) 522 522 ! 523 !-- Limit dissipation rate according to Seifert, Nuijens and 523 !-- Limit dissipation rate according to Seifert, Nuijens and 524 524 !-- Stevens (2010) 525 525 dissipation = MIN( 0.06_wp, diss(k,j,i) ) … … 531 531 dissipation**( 1.0_wp / 6.0_wp ) 532 532 ! 533 !-- The factor of 1.0E4 is needed to convert the dissipation 533 !-- The factor of 1.0E4 is needed to convert the dissipation 534 534 !-- rate from m2 s-3 to cm2 s-3. 535 535 k_au = k_au * ( 1.0_wp + & 536 536 dissipation * 1.0E4_wp * & 537 537 ( re_lambda * 1.0E-3_wp )**0.25_wp * & 538 ( alpha_cc * EXP( -1.0_wp * ( ( rc - & 538 ( alpha_cc * EXP( -1.0_wp * ( ( rc - & 539 539 r_cc ) / & 540 540 sigma_cc )**2 & … … 552 552 553 553 qr(k,j,i) = qr(k,j,i) + autocon * dt_micro 554 qc(k,j,i) = qc(k,j,i) - autocon * dt_micro 554 qc(k,j,i) = qc(k,j,i) - autocon * dt_micro 555 555 nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro 556 556 … … 580 580 581 581 USE indices, & 582 ONLY: nxl , nxr, nyn, nys, nzb_s_inner, nzt582 ONLY: nxlg, nxrg, nyng, nysg, nzb_s_inner, nzt 583 583 584 584 USE kinds … … 593 593 REAL(wp) :: dqdt_precip !< 594 594 595 DO i = nxl , nxr596 DO j = nys , nyn595 DO i = nxlg, nxrg 596 DO j = nysg, nyng 597 597 DO k = nzb_s_inner(j,i)+1, nzt 598 598 … … 640 640 641 641 USE indices, & 642 ONLY: nxl , nxr, nys, nyn, nzb_s_inner, nzt642 ONLY: nxlg, nxrg, nysg, nyng, nzb_s_inner, nzt 643 643 644 644 USE kinds … … 657 657 CALL cpu_log( log_point_s(56), 'accretion', 'start' ) 658 658 659 DO i = nxl , nxr660 DO j = nys , nyn659 DO i = nxlg, nxrg 660 DO j = nysg, nyng 661 661 DO k = nzb_s_inner(j,i)+1, nzt 662 662 … … 664 664 ! 665 665 !-- Intern time scale of coagulation (Seifert and Beheng, 2006): 666 tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ) 667 ! 668 !-- Universal function for accretion process (Seifert and 666 tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ) 667 ! 668 !-- Universal function for accretion process (Seifert and 669 669 !-- Beheng, 2001): 670 670 phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4 671 671 ! 672 !-- Parameterized turbulence effects on autoconversion (Seifert, 673 !-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to 672 !-- Parameterized turbulence effects on autoconversion (Seifert, 673 !-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to 674 674 !-- convert the dissipation rate (diss) from m2 s-3 to cm2 s-3. 675 675 IF ( collision_turbulence ) THEN … … 679 679 ) 680 680 ELSE 681 k_cr = k_cr0 681 k_cr = k_cr0 682 682 ENDIF 683 683 ! … … 687 687 accr = MIN( accr, qc(k,j,i) / dt_micro ) 688 688 689 qr(k,j,i) = qr(k,j,i) + accr * dt_micro 689 qr(k,j,i) = qr(k,j,i) + accr * dt_micro 690 690 qc(k,j,i) = qc(k,j,i) - accr * dt_micro 691 691 … … 721 721 722 722 USE indices, & 723 ONLY: nxl , nxr, nys, nyn, nzb_s_inner, nzt723 ONLY: nxlg, nxrg, nysg, nyng, nzb_s_inner, nzt 724 724 725 725 USE kinds 726 726 727 727 IMPLICIT NONE 728 728 … … 738 738 CALL cpu_log( log_point_s(57), 'selfcollection', 'start' ) 739 739 740 DO i = nxl , nxr741 DO j = nys , nyn740 DO i = nxlg, nxrg 741 DO j = nysg, nyng 742 742 DO k = nzb_s_inner(j,i)+1, nzt 743 743 IF ( qr(k,j,i) > eps_sb ) THEN … … 762 762 nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro 763 763 764 ENDIF 764 ENDIF 765 765 ENDDO 766 766 ENDDO … … 775 775 ! Description: 776 776 ! ------------ 777 !> Evaporation of precipitable water. Condensation is neglected for 777 !> Evaporation of precipitable water. Condensation is neglected for 778 778 !> precipitable water. 779 779 !------------------------------------------------------------------------------! … … 793 793 794 794 USE indices, & 795 ONLY: nxl , nxr, nys, nyn, nzb_s_inner, nzt795 ONLY: nxlg, nxrg, nysg, nyng, nzb_s_inner, nzt 796 796 797 797 USE kinds … … 823 823 CALL cpu_log( log_point_s(58), 'evaporation', 'start' ) 824 824 825 DO i = nxl , nxr826 DO j = nys , nyn825 DO i = nxlg, nxrg 826 DO j = nysg, nyng 827 827 DO k = nzb_s_inner(j,i)+1, nzt 828 828 IF ( qr(k,j,i) > eps_sb ) THEN … … 850 850 !-- Actual temperature: 851 851 temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) ) 852 852 853 853 g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * & 854 854 l_v / ( thermal_conductivity_l * temp ) & … … 862 862 dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp ) 863 863 ! 864 !-- Compute ventilation factor and intercept parameter 864 !-- Compute ventilation factor and intercept parameter 865 865 !-- (Seifert and Beheng, 2006; Seifert, 2008): 866 866 IF ( ventilation_effect ) THEN 867 867 ! 868 !-- Shape parameter of gamma distribution (Milbrandt and Yau, 868 !-- Shape parameter of gamma distribution (Milbrandt and Yau, 869 869 !-- 2005; Stevens and Seifert, 2008): 870 870 mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * & … … 877 877 878 878 mu_r_2 = mu_r + 2.0_wp 879 mu_r_5d2 = mu_r + 2.5_wp 879 mu_r_5d2 = mu_r + 2.5_wp 880 880 881 881 f_vent = a_vent * gamm( mu_r_2 ) * & … … 893 893 ( lambda_r / ( 3.0_wp * c_term + lambda_r ) & 894 894 )**mu_r_5d2 - & 895 0.0390625_wp * ( b_term / a_term )**4 * & 895 0.0390625_wp * ( b_term / a_term )**4 * & 896 896 ( lambda_r / ( 4.0_wp * c_term + lambda_r ) & 897 897 )**mu_r_5d2 & … … 899 899 900 900 nr_0 = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) / & 901 gamm( mu_r + 1.0_wp ) 901 gamm( mu_r + 1.0_wp ) 902 902 ELSE 903 903 f_vent = 1.0_wp … … 905 905 ENDIF 906 906 ! 907 !-- Evaporation rate of rain water content (Seifert and 907 !-- Evaporation rate of rain water content (Seifert and 908 908 !-- Beheng, 2006): 909 909 evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / & … … 917 917 918 918 ENDIF 919 ENDIF 919 ENDIF 920 920 921 921 ENDDO … … 948 948 949 949 USE indices, & 950 ONLY: nxl , nxr, nys, nyn, nzb, nzb_s_inner, nzt950 ONLY: nxlg, nxrg, nysg, nyng, nzb, nzb_s_inner, nzt 951 951 952 952 USE kinds … … 954 954 USE statistics, & 955 955 ONLY: weight_substep 956 956 957 957 958 958 IMPLICIT NONE … … 968 968 sed_qc(nzt+1) = 0.0_wp 969 969 970 DO i = nxl , nxr971 DO j = nys , nyn970 DO i = nxlg, nxrg 971 DO j = nysg, nyng 972 972 DO k = nzt, nzb_s_inner(j,i)+1, -1 973 973 … … 1029 1029 1030 1030 USE indices, & 1031 ONLY: nxl , nxr, nys, nyn, nzb, nzb_s_inner, nzt1031 ONLY: nxlg, nxrg, nysg, nyng, nzb, nzb_s_inner, nzt 1032 1032 1033 1033 USE kinds … … 1035 1035 USE statistics, & 1036 1036 ONLY: weight_substep 1037 1037 1038 1038 IMPLICIT NONE 1039 1039 … … 1065 1065 1066 1066 ! 1067 !-- Compute velocities 1068 DO i = nxl , nxr1069 DO j = nys , nyn1067 !-- Compute velocities 1068 DO i = nxlg, nxrg 1069 DO j = nysg, nyng 1070 1070 DO k = nzb_s_inner(j,i)+1, nzt 1071 1071 IF ( qr(k,j,i) > eps_sb ) THEN … … 1119 1119 2.0_wp * w_qr(k) + w_qr(k+1) ) * & 1120 1120 dt_micro * ddzu(k) 1121 ENDDO 1121 ENDDO 1122 1122 ! 1123 1123 !-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977): … … 1155 1155 DO k = nzt, nzb_s_inner(j,i)+1, -1 1156 1156 ! 1157 !-- Sum up all rain drop number densities which contribute to the flux 1157 !-- Sum up all rain drop number densities which contribute to the flux 1158 1158 !-- through k-1/2 1159 1159 flux = 0.0_wp … … 1170 1170 ENDDO 1171 1171 ! 1172 !-- It is not allowed to sediment more rain drop number density than 1172 !-- It is not allowed to sediment more rain drop number density than 1173 1173 !-- available 1174 1174 flux = MIN( flux, & … … 1181 1181 ddzu(k+1) / hyrho(k) * dt_micro 1182 1182 ! 1183 !-- Sum up all rain water content which contributes to the flux 1183 !-- Sum up all rain water content which contributes to the flux 1184 1184 !-- through k-1/2 1185 1185 flux = 0.0_wp … … 1199 1199 ENDDO 1200 1200 ! 1201 !-- It is not allowed to sediment more rain water content than 1201 !-- It is not allowed to sediment more rain water content than 1202 1202 !-- available 1203 1203 flux = MIN( flux, & … … 1211 1211 ddzu(k+1) / hyrho(k) * dt_micro 1212 1212 q(k,j,i) = q(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * & 1213 ddzu(k+1) / hyrho(k) * dt_micro 1213 ddzu(k+1) / hyrho(k) * dt_micro 1214 1214 pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) * & 1215 1215 ddzu(k+1) / hyrho(k) * l_d_cp * & … … 1330 1330 pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp 1331 1331 t_d_pt(k) = 1.0_wp / pt_d_t(k) 1332 hyrho(k) = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) ) 1332 hyrho(k) = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) ) 1333 1333 ENDDO 1334 1334 ! … … 1338 1338 1339 1339 ! 1340 !-- Compute length of time step 1340 !-- Compute length of time step 1341 1341 IF ( call_microphysics_at_all_substeps ) THEN 1342 1342 dt_micro = dt_3d * weight_pres(intermediate_timestep_count) … … 1395 1395 ! Description: 1396 1396 ! ------------ 1397 !> Adjust number of raindrops to avoid nonlinear effects in 1397 !> Adjust number of raindrops to avoid nonlinear effects in 1398 1398 !> sedimentation and evaporation of rain drops due to too small or 1399 1399 !> too big weights of rain drops (Stevens and Seifert, 2008). … … 1424 1424 ELSE 1425 1425 ! 1426 !-- Adjust number of raindrops to avoid nonlinear effects in 1426 !-- Adjust number of raindrops to avoid nonlinear effects in 1427 1427 !-- sedimentation and evaporation of rain drops due to too small or 1428 1428 !-- too big weights of rain drops (Stevens and Seifert, 2008). … … 1470 1470 INTEGER(iwp) :: k !< 1471 1471 1472 REAL(wp) :: alpha_cc !< 1472 REAL(wp) :: alpha_cc !< 1473 1473 REAL(wp) :: autocon !< 1474 1474 REAL(wp) :: dissipation !< … … 1494 1494 tau_cloud = 1.0_wp - qc_1d(k) / ( qr_1d(k) + qc_1d(k) ) 1495 1495 ! 1496 !-- Universal function for autoconversion process 1496 !-- Universal function for autoconversion process 1497 1497 !-- (Seifert and Beheng, 2006): 1498 1498 phi_au = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**3 … … 1505 1505 xc = hyrho(k) * qc_1d(k) / nc_1d(k) 1506 1506 ! 1507 !-- Parameterized turbulence effects on autoconversion (Seifert, 1507 !-- Parameterized turbulence effects on autoconversion (Seifert, 1508 1508 !-- Nuijens and Stevens, 2010) 1509 1509 IF ( collision_turbulence ) THEN … … 1519 1519 l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp ) 1520 1520 ! 1521 !-- Limit dissipation rate according to Seifert, Nuijens and 1521 !-- Limit dissipation rate according to Seifert, Nuijens and 1522 1522 !-- Stevens (2010) 1523 1523 dissipation = MIN( 0.06_wp, diss(k,j,i) ) … … 1549 1549 1550 1550 qr_1d(k) = qr_1d(k) + autocon * dt_micro 1551 qc_1d(k) = qc_1d(k) - autocon * dt_micro 1551 qc_1d(k) = qc_1d(k) - autocon * dt_micro 1552 1552 nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro 1553 1553 … … 1642 1642 ! 1643 1643 !-- Intern time scale of coagulation (Seifert and Beheng, 2006): 1644 tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) ) 1645 ! 1646 !-- Universal function for accretion process 1644 tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) ) 1645 ! 1646 !-- Universal function for accretion process 1647 1647 !-- (Seifert and Beheng, 2001): 1648 1648 phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4 1649 1649 ! 1650 !-- Parameterized turbulence effects on autoconversion (Seifert, 1651 !-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to 1650 !-- Parameterized turbulence effects on autoconversion (Seifert, 1651 !-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to 1652 1652 !-- convert the dissipation rate (diss) from m2 s-3 to cm2 s-3. 1653 1653 IF ( collision_turbulence ) THEN … … 1657 1657 ) 1658 1658 ELSE 1659 k_cr = k_cr0 1659 k_cr = k_cr0 1660 1660 ENDIF 1661 1661 ! … … 1664 1664 accr = MIN( accr, qc_1d(k) / dt_micro ) 1665 1665 1666 qr_1d(k) = qr_1d(k) + accr * dt_micro 1666 qr_1d(k) = qr_1d(k) + accr * dt_micro 1667 1667 qc_1d(k) = qc_1d(k) - accr * dt_micro 1668 1668 … … 1691 1691 1692 1692 USE kinds 1693 1693 1694 1694 IMPLICIT NONE 1695 1695 … … 1723 1723 nr_1d(k) = nr_1d(k) + selfcoll * dt_micro 1724 1724 1725 ENDIF 1725 ENDIF 1726 1726 ENDDO 1727 1727 … … 1732 1732 ! Description: 1733 1733 ! ------------ 1734 !> Evaporation of precipitable water. Condensation is neglected for 1734 !> Evaporation of precipitable water. Condensation is neglected for 1735 1735 !> precipitable water. Call for grid point i,j 1736 1736 !------------------------------------------------------------------------------! … … 1799 1799 !-- Actual temperature: 1800 1800 temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) ) 1801 1801 1802 1802 g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v / & 1803 1803 ( thermal_conductivity_l * temp ) + & … … 1811 1811 dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp ) 1812 1812 ! 1813 !-- Compute ventilation factor and intercept parameter 1813 !-- Compute ventilation factor and intercept parameter 1814 1814 !-- (Seifert and Beheng, 2006; Seifert, 2008): 1815 1815 IF ( ventilation_effect ) THEN … … 1825 1825 1826 1826 mu_r_2 = mu_r + 2.0_wp 1827 mu_r_5d2 = mu_r + 2.5_wp 1828 1829 f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) + & 1827 mu_r_5d2 = mu_r + 2.5_wp 1828 1829 f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) + & 1830 1830 b_vent * schmidt_p_1d3 * & 1831 1831 SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) * & … … 1841 1841 ( lambda_r / ( 3.0_wp * c_term + lambda_r ) & 1842 1842 )**mu_r_5d2 - & 1843 0.0390625_wp * ( b_term / a_term )**4 * & 1843 0.0390625_wp * ( b_term / a_term )**4 * & 1844 1844 ( lambda_r / ( 4.0_wp * c_term + lambda_r ) & 1845 1845 )**mu_r_5d2 & … … 1847 1847 1848 1848 nr_0 = nr_1d(k) * lambda_r**( mu_r + 1.0_wp ) / & 1849 gamm( mu_r + 1.0_wp ) 1849 gamm( mu_r + 1.0_wp ) 1850 1850 ELSE 1851 1851 f_vent = 1.0_wp … … 1863 1863 1864 1864 ENDIF 1865 ENDIF 1865 ENDIF 1866 1866 1867 1867 ENDDO … … 1873 1873 ! Description: 1874 1874 ! ------------ 1875 !> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR). 1875 !> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR). 1876 1876 !> Call for grid point i,j 1877 1877 !------------------------------------------------------------------------------! … … 1891 1891 1892 1892 USE kinds 1893 1893 1894 1894 USE statistics, & 1895 1895 ONLY: weight_substep … … 1919 1919 q_1d(k) = q_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / & 1920 1920 hyrho(k) * dt_micro 1921 qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / & 1921 qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / & 1922 1922 hyrho(k) * dt_micro 1923 1923 pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / & … … 1962 1962 USE statistics, & 1963 1963 ONLY: weight_substep 1964 1964 1965 1965 IMPLICIT NONE 1966 1966 … … 1990 1990 1991 1991 ! 1992 !-- Compute velocities 1992 !-- Compute velocities 1993 1993 DO k = nzb_s_inner(j,i)+1, nzt 1994 1994 IF ( qr_1d(k) > eps_sb ) THEN … … 2035 2035 c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) * & 2036 2036 dt_micro * ddzu(k) 2037 ENDDO 2037 ENDDO 2038 2038 ! 2039 2039 !-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977): … … 2071 2071 DO k = nzt, nzb_s_inner(j,i)+1, -1 2072 2072 ! 2073 !-- Sum up all rain drop number densities which contribute to the flux 2073 !-- Sum up all rain drop number densities which contribute to the flux 2074 2074 !-- through k-1/2 2075 2075 flux = 0.0_wp … … 2086 2086 ENDDO 2087 2087 ! 2088 !-- It is not allowed to sediment more rain drop number density than 2088 !-- It is not allowed to sediment more rain drop number density than 2089 2089 !-- available 2090 2090 flux = MIN( flux, & … … 2095 2095 hyrho(k) * dt_micro 2096 2096 ! 2097 !-- Sum up all rain water content which contributes to the flux 2097 !-- Sum up all rain water content which contributes to the flux 2098 2098 !-- through k-1/2 2099 2099 flux = 0.0_wp … … 2122 2122 hyrho(k) * dt_micro 2123 2123 q_1d(k) = q_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / & 2124 hyrho(k) * dt_micro 2124 hyrho(k) * dt_micro 2125 2125 pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / & 2126 2126 hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro … … 2185 2185 ! ------------ 2186 2186 !> This function computes the gamma function (Press et al., 1992). 2187 !> The gamma function is needed for the calculation of the evaporation 2187 !> The gamma function is needed for the calculation of the evaporation 2188 2188 !> of rain drops. 2189 2189 !------------------------------------------------------------------------------! 2190 FUNCTION gamm( xx ) 2190 FUNCTION gamm( xx ) 2191 2191 2192 2192 USE kinds 2193 2193 2194 IMPLICIT NONE 2194 IMPLICIT NONE 2195 2195 2196 2196 INTEGER(iwp) :: j !< … … 2218 2218 ser = 1.000000000190015_wp 2219 2219 2220 DO j = 1, 6 2221 y_gamm = y_gamm + 1.0_wp 2222 ser = ser + cof( j ) / y_gamm 2220 DO j = 1, 6 2221 y_gamm = y_gamm + 1.0_wp 2222 ser = ser + cof( j ) / y_gamm 2223 2223 ENDDO 2224 2224 2225 ! 2226 !-- Until this point the algorithm computes the logarithm of the gamma 2227 !-- function. Hence, the exponential function is used. 2228 ! gamm = EXP( tmp + LOG( stp * ser / x_gamm ) ) 2229 gamm = EXP( tmp ) * stp * ser / x_gamm 2230 2231 RETURN 2232 2233 END FUNCTION gamm 2225 ! 2226 !-- Until this point the algorithm computes the logarithm of the gamma 2227 !-- function. Hence, the exponential function is used. 2228 ! gamm = EXP( tmp + LOG( stp * ser / x_gamm ) ) 2229 gamm = EXP( tmp ) * stp * ser / x_gamm 2230 2231 RETURN 2232 2233 END FUNCTION gamm 2234 2234 2235 2235 END MODULE microphysics_mod
Note: See TracChangeset
for help on using the changeset viewer.
