Changeset 2155 for palm/trunk
- Timestamp:
- Feb 21, 2017 9:57:40 AM (8 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 2 edited
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 -
palm/trunk/SOURCE/prognostic_equations.f90
r2119 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 ! 2118 2017-01-17 16:38:49Z raasch 29 29 ! OpenACC version of subroutine removed 30 ! 30 ! 31 31 ! 2031 2016-10-21 15:11:58Z knoop 32 32 ! renamed variable rho to rho_ocean 33 ! 33 ! 34 34 ! 2011 2016-09-19 17:29:57Z kanani 35 35 ! Flag urban_surface is now defined in module control_parameters. 36 ! 36 ! 37 37 ! 2007 2016-08-24 15:47:17Z kanani 38 38 ! Added pt tendency calculation based on energy balance at urban surfaces 39 39 ! (new urban surface model) 40 ! 40 ! 41 41 ! 2000 2016-08-20 18:09:15Z knoop 42 42 ! Forced header and separation lines into 80 columns 43 ! 43 ! 44 44 ! 1976 2016-07-27 13:28:04Z maronga 45 45 ! Simplied calls to radiation model 46 ! 46 ! 47 47 ! 1960 2016-07-12 16:34:24Z suehring 48 48 ! Separate humidity and passive scalar 49 ! 49 ! 50 50 ! 1914 2016-05-26 14:44:07Z witha 51 51 ! Added calls for wind turbine model … … 53 53 ! 1873 2016-04-18 14:50:06Z maronga 54 54 ! Module renamed (removed _mod) 55 ! 55 ! 56 56 ! 1850 2016-04-08 13:29:27Z maronga 57 57 ! Module renamed 58 ! 58 ! 59 59 ! 1826 2016-04-07 12:01:39Z maronga 60 60 ! Renamed canopy model calls. 61 ! 61 ! 62 62 ! 1822 2016-04-07 07:49:42Z hoffmann 63 63 ! Kessler microphysics scheme moved to microphysics. 64 64 ! 65 65 ! 1757 2016-02-22 15:49:32Z maronga 66 ! 66 ! 67 67 ! 1691 2015-10-26 16:17:44Z maronga 68 68 ! Added optional model spin-up without radiation / land surface model calls. 69 69 ! Formatting corrections. 70 ! 70 ! 71 71 ! 1682 2015-10-07 23:56:08Z knoop 72 ! Code annotations made doxygen readable 73 ! 72 ! Code annotations made doxygen readable 73 ! 74 74 ! 1585 2015-04-30 07:05:52Z maronga 75 75 ! Added call for temperature tendency calculation due to radiative flux divergence 76 ! 76 ! 77 77 ! 1517 2015-01-07 19:12:25Z hoffmann 78 78 ! advec_s_bc_mod addded, since advec_s_bc is now a module … … 80 80 ! 1496 2014-12-02 17:25:50Z maronga 81 81 ! Renamed "radiation" -> "cloud_top_radiation" 82 ! 82 ! 83 83 ! 1484 2014-10-21 10:53:05Z kanani 84 84 ! Changes due to new module structure of the plant canopy model: … … 86 86 ! Removed double-listing of use_upstream_for_tke in ONLY-list of module 87 87 ! control_parameters 88 ! 88 ! 89 89 ! 1409 2014-05-23 12:11:32Z suehring 90 ! Bugfix: i_omp_start changed for advec_u_ws at left inflow and outflow boundary. 90 ! Bugfix: i_omp_start changed for advec_u_ws at left inflow and outflow boundary. 91 91 ! This ensures that left-hand side fluxes are also calculated for nxl in that 92 ! case, even though the solution at nxl is overwritten in boundary_conds() 93 ! 92 ! case, even though the solution at nxl is overwritten in boundary_conds() 93 ! 94 94 ! 1398 2014-05-07 11:15:00Z heinze 95 95 ! Rayleigh-damping for horizontal velocity components changed: instead of damping 96 ! against ug and vg, damping against u_init and v_init is used to allow for a 96 ! against ug and vg, damping against u_init and v_init is used to allow for a 97 97 ! homogenized treatment in case of nudging 98 ! 98 ! 99 99 ! 1380 2014-04-28 12:40:45Z heinze 100 ! Change order of calls for scalar prognostic quantities: 101 ! ls_advec -> nudging -> subsidence since initial profiles 102 ! 100 ! Change order of calls for scalar prognostic quantities: 101 ! ls_advec -> nudging -> subsidence since initial profiles 102 ! 103 103 ! 1374 2014-04-25 12:55:07Z raasch 104 104 ! missing variables added to ONLY lists 105 ! 105 ! 106 106 ! 1365 2014-04-22 15:03:56Z boeske 107 ! Calls of ls_advec for large scale advection added, 107 ! Calls of ls_advec for large scale advection added, 108 108 ! subroutine subsidence is only called if use_subsidence_tendencies = .F., 109 109 ! new argument ls_index added to the calls of subsidence 110 110 ! +ls_index 111 ! 111 ! 112 112 ! 1361 2014-04-16 15:17:48Z hoffmann 113 113 ! Two-moment microphysics moved to the start of prognostic equations. This makes … … 117 117 ! 118 118 ! Two-moment cloud physics added for vector and accelerator optimization. 119 ! 119 ! 120 120 ! 1353 2014-04-08 15:21:23Z heinze 121 121 ! REAL constants provided with KIND-attribute 122 ! 122 ! 123 123 ! 1337 2014-03-25 15:11:48Z heinze 124 124 ! Bugfix: REAL constants provided with KIND-attribute 125 ! 125 ! 126 126 ! 1332 2014-03-25 11:59:43Z suehring 127 ! Bugfix: call advec_ws or advec_pw for TKE only if NOT use_upstream_for_tke 128 ! 127 ! Bugfix: call advec_ws or advec_pw for TKE only if NOT use_upstream_for_tke 128 ! 129 129 ! 1330 2014-03-24 17:29:32Z suehring 130 ! In case of SGS-particle velocity advection of TKE is also allowed with 130 ! In case of SGS-particle velocity advection of TKE is also allowed with 131 131 ! dissipative 5th-order scheme. 132 132 ! … … 161 161 ! 162 162 ! 1115 2013-03-26 18:16:16Z hoffmann 163 ! optimized cloud physics: calculation of microphysical tendencies transfered 163 ! optimized cloud physics: calculation of microphysical tendencies transfered 164 164 ! to microphysics.f90; qr and nr are only calculated if precipitation is required 165 165 ! … … 211 211 MODULE prognostic_equations_mod 212 212 213 213 214 214 215 215 USE arrays_3d, & … … 225 225 te_m, tnr_m, tpt_m, tq_m, tqr_m, ts_m, tsa_m, tswst, tu_m, tv_m,& 226 226 tw_m, u, ug, u_init, u_p, v, vg, vpt, v_init, v_p, w, w_p 227 227 228 228 USE control_parameters, & 229 229 ONLY: call_microphysics_at_all_substeps, cloud_physics, & … … 249 249 250 250 USE indices, & 251 ONLY: nxl, nxl u, nxr, nyn, nys, nysv, nzb_s_inner, nzb_u_inner,&252 nzb_ v_inner, nzb_w_inner, nzt251 ONLY: nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv, & 252 nzb_s_inner, nzb_u_inner, nzb_v_inner, nzb_w_inner, nzt 253 253 254 254 USE advec_ws, & … … 287 287 USE calc_radiation_mod, & 288 288 ONLY: calc_radiation 289 289 290 290 USE coriolis_mod, & 291 291 ONLY: coriolis … … 370 370 !> (Optimized to avoid cache missings, i.e. for Power4/5-architectures.) 371 371 !------------------------------------------------------------------------------! 372 372 373 373 SUBROUTINE prognostic_equations_cache 374 374 … … 382 382 INTEGER(iwp) :: omp_get_thread_num !< 383 383 INTEGER(iwp) :: tn = 0 !< 384 384 385 385 LOGICAL :: loop_start !< 386 386 … … 394 394 !$OMP PARALLEL private (i,i_omp_start,j,k,loop_start,tn) 395 395 396 !$ tn = omp_get_thread_num() 396 !$ tn = omp_get_thread_num() 397 397 loop_start = .TRUE. 398 398 !$OMP DO 399 400 ! 401 !-- If required, calculate cloud microphysics 402 IF ( cloud_physics .AND. .NOT. microphysics_sat_adjust .AND. & 403 ( intermediate_timestep_count == 1 .OR. & 404 call_microphysics_at_all_substeps ) & 405 ) THEN 406 DO i = nxlg, nxrg 407 DO j = nysg, nyng 408 CALL microphysics_control( i, j ) 409 END DO 410 END DO 411 ENDIF 412 399 413 DO i = nxl, nxr 400 414 … … 404 418 IF ( loop_start ) THEN 405 419 loop_start = .FALSE. 406 i_omp_start = i 420 i_omp_start = i 407 421 ENDIF 408 422 409 423 DO j = nys, nyn 410 !411 !-- If required, calculate cloud microphysics412 IF ( cloud_physics .AND. .NOT. microphysics_sat_adjust .AND. &413 ( intermediate_timestep_count == 1 .OR. &414 call_microphysics_at_all_substeps ) &415 ) THEN416 CALL microphysics_control( i, j )417 ENDIF418 424 ! 419 425 !-- Tendency terms for u-velocity component … … 424 430 IF ( ws_scheme_mom ) THEN 425 431 CALL advec_u_ws( i, j, i_omp_start, tn ) 426 ELSE 432 ELSE 427 433 CALL advec_u_pw( i, j ) 428 ENDIF 434 ENDIF 429 435 ELSE 430 436 CALL advec_u_up( i, j ) … … 490 496 IF ( ws_scheme_mom ) THEN 491 497 CALL advec_v_ws( i, j, i_omp_start, tn ) 492 ELSE 498 ELSE 493 499 CALL advec_v_pw( i, j ) 494 500 ENDIF … … 501 507 ! 502 508 !-- Drag by plant canopy 503 IF ( plant_canopy ) CALL pcm_tendency( i, j, 2 ) 509 IF ( plant_canopy ) CALL pcm_tendency( i, j, 2 ) 504 510 505 511 ! … … 551 557 IF ( ws_scheme_mom ) THEN 552 558 CALL advec_w_ws( i, j, i_omp_start, tn ) 553 ELSE 559 ELSE 554 560 CALL advec_w_pw( i, j ) 555 561 END IF … … 646 652 IF ( large_scale_forcing ) THEN 647 653 CALL ls_advec( i, j, simulated_time, 'pt' ) 648 ENDIF 654 ENDIF 649 655 650 656 ! 651 657 !-- Nudging 652 IF ( nudging ) CALL nudge( i, j, simulated_time, 'pt' ) 658 IF ( nudging ) CALL nudge( i, j, simulated_time, 'pt' ) 653 659 654 660 ! … … 708 714 CALL advec_s_ws( i, j, sa, 'sa', flux_s_sa, & 709 715 diss_s_sa, flux_l_sa, diss_l_sa, i_omp_start, tn ) 710 ELSE 716 ELSE 711 717 CALL advec_s_pw( i, j, sa ) 712 718 ENDIF … … 760 766 THEN 761 767 IF ( ws_scheme_sca ) THEN 762 CALL advec_s_ws( i, j, q, 'q', flux_s_q, & 768 CALL advec_s_ws( i, j, q, 'q', flux_s_q, & 763 769 diss_s_q, flux_l_q, diss_l_q, i_omp_start, tn ) 764 ELSE 770 ELSE 765 771 CALL advec_s_pw( i, j, q ) 766 772 ENDIF … … 782 788 ! 783 789 !-- Nudging 784 IF ( nudging ) CALL nudge( i, j, simulated_time, 'q' ) 790 IF ( nudging ) CALL nudge( i, j, simulated_time, 'q' ) 785 791 786 792 ! … … 820 826 821 827 ! 822 !-- If required, calculate prognostic equations for rain water content 828 !-- If required, calculate prognostic equations for rain water content 823 829 !-- and rain drop concentration 824 830 IF ( cloud_physics .AND. microphysics_seifert ) THEN … … 829 835 THEN 830 836 IF ( ws_scheme_sca ) THEN 831 CALL advec_s_ws( i, j, qr, 'qr', flux_s_qr, & 837 CALL advec_s_ws( i, j, qr, 'qr', flux_s_qr, & 832 838 diss_s_qr, flux_l_qr, diss_l_qr, & 833 839 i_omp_start, tn ) 834 ELSE 840 ELSE 835 841 CALL advec_s_pw( i, j, qr ) 836 842 ENDIF … … 869 875 IF ( timestep_scheme(1:5) == 'runge' ) THEN 870 876 IF ( ws_scheme_sca ) THEN 871 CALL advec_s_ws( i, j, nr, 'nr', flux_s_nr, & 877 CALL advec_s_ws( i, j, nr, 'nr', flux_s_nr, & 872 878 diss_s_nr, flux_l_nr, diss_l_nr, & 873 879 i_omp_start, tn ) 874 ELSE 880 ELSE 875 881 CALL advec_s_pw( i, j, nr ) 876 882 ENDIF … … 907 913 908 914 ENDIF 909 915 910 916 ! 911 917 !-- If required, compute prognostic equation for scalar … … 917 923 THEN 918 924 IF ( ws_scheme_sca ) THEN 919 CALL advec_s_ws( i, j, s, 's', flux_s_s, & 925 CALL advec_s_ws( i, j, s, 's', flux_s_s, & 920 926 diss_s_s, flux_l_s, diss_l_s, i_omp_start, tn ) 921 ELSE 927 ELSE 922 928 CALL advec_s_pw( i, j, s ) 923 929 ENDIF … … 939 945 ! 940 946 !-- Nudging, still need to be extended for scalars 941 ! IF ( nudging ) CALL nudge( i, j, simulated_time, 's' ) 947 ! IF ( nudging ) CALL nudge( i, j, simulated_time, 's' ) 942 948 943 949 ! 944 950 !-- If required compute influence of large-scale subsidence/ascent. 945 !-- Note, the last argument is of no meaning in this case, as it is 946 !-- only used in conjunction with large_scale_forcing, which is to 951 !-- Note, the last argument is of no meaning in this case, as it is 952 !-- only used in conjunction with large_scale_forcing, which is to 947 953 !-- date not implemented for scalars. 948 954 IF ( large_scale_subsidence .AND. & … … 980 986 ENDIF 981 987 982 ENDIF 983 ! 984 !-- If required, compute prognostic equation for turbulent kinetic 988 ENDIF 989 ! 990 !-- If required, compute prognostic equation for turbulent kinetic 985 991 !-- energy (TKE) 986 992 IF ( .NOT. constant_diffusion ) THEN … … 990 996 tend(:,j,i) = 0.0_wp 991 997 IF ( timestep_scheme(1:5) == 'runge' & 992 .AND. .NOT. use_upstream_for_tke ) THEN 998 .AND. .NOT. use_upstream_for_tke ) THEN 993 999 IF ( ws_scheme_sca ) THEN 994 1000 CALL advec_s_ws( i, j, e, 'e', flux_s_e, diss_s_e, & … … 1013 1019 ! 1014 1020 !-- Additional sink term for flows through plant canopies 1015 IF ( plant_canopy ) CALL pcm_tendency( i, j, 6 ) 1021 IF ( plant_canopy ) CALL pcm_tendency( i, j, 6 ) 1016 1022 1017 1023 CALL user_actions( i, j, 'e-tendency' ) … … 1061 1067 !> Version for vector machines 1062 1068 !------------------------------------------------------------------------------! 1063 1069 1064 1070 SUBROUTINE prognostic_equations_vector 1065 1071 … … 1176 1182 IF ( ws_scheme_mom ) THEN 1177 1183 CALL advec_v_ws 1178 ELSE 1184 ELSE 1179 1185 CALL advec_v_pw 1180 1186 END IF … … 1389 1395 ! 1390 1396 !-- Nudging 1391 IF ( nudging ) CALL nudge( simulated_time, 'pt' ) 1397 IF ( nudging ) CALL nudge( simulated_time, 'pt' ) 1392 1398 1393 1399 ! … … 1485 1491 1486 1492 CALL diffusion_s( sa, saswsb, saswst, wall_salinityflux ) 1487 1493 1488 1494 CALL user_actions( 'sa-tendency' ) 1489 1495 … … 1573 1579 1574 1580 CALL diffusion_s( q, qsws, qswst, wall_qflux ) 1575 1581 1576 1582 ! 1577 1583 !-- Sink or source of humidity due to canopy elements … … 1586 1592 ! 1587 1593 !-- Nudging 1588 IF ( nudging ) CALL nudge( simulated_time, 'q' ) 1594 IF ( nudging ) CALL nudge( simulated_time, 'q' ) 1589 1595 1590 1596 ! … … 1637 1643 1638 1644 ! 1639 !-- If required, calculate prognostic equations for rain water content 1645 !-- If required, calculate prognostic equations for rain water content 1640 1646 !-- and rain drop concentration 1641 1647 IF ( cloud_physics .AND. microphysics_seifert ) THEN … … 1675 1681 CALL diffusion_s( qr, qrsws, qrswst, wall_qrflux ) 1676 1682 1677 CALL user_actions( 'qr-tendency' )1678 1679 1683 ! 1680 1684 !-- Prognostic equation for rain water content … … 1828 1832 1829 1833 CALL diffusion_s( s, ssws, sswst, wall_sflux ) 1830 1834 1831 1835 ! 1832 1836 !-- Sink or source of humidity due to canopy elements … … 1841 1845 ! 1842 1846 !-- Nudging. Not implemented for scalars so far. 1843 ! IF ( nudging ) CALL nudge( simulated_time, 'q' ) 1847 ! IF ( nudging ) CALL nudge( simulated_time, 'q' ) 1844 1848 1845 1849 ! … … 1895 1899 ENDIF 1896 1900 ! 1897 !-- If required, compute prognostic equation for turbulent kinetic 1901 !-- If required, compute prognostic equation for turbulent kinetic 1898 1902 !-- energy (TKE) 1899 1903 IF ( .NOT. constant_diffusion ) THEN … … 1927 1931 IF ( ws_scheme_sca ) THEN 1928 1932 CALL advec_s_ws( e, 'e' ) 1929 ELSE 1933 ELSE 1930 1934 CALL advec_s_pw( e ) 1931 1935 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.