Changeset 1822 for palm/trunk/SOURCE/microphysics.f90
- Timestamp:
- Apr 7, 2016 7:49:42 AM (8 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/microphysics.f90 (modified) (45 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/microphysics.f90
r1818 r1822 19 19 ! Current revisions: 20 20 ! ------------------ 21 ! 21 ! Unused variables removed. 22 ! 23 ! Kessler scheme integrated. 22 24 ! 23 25 ! Former revisions: … … 109 111 END INTERFACE autoconversion 110 112 113 INTERFACE autoconversion_kessler 114 MODULE PROCEDURE autoconversion_kessler 115 MODULE PROCEDURE autoconversion_kessler_ij 116 END INTERFACE autoconversion_kessler 117 111 118 INTERFACE accretion 112 119 MODULE PROCEDURE accretion … … 151 158 152 159 USE arrays_3d, & 153 ONLY: hyp, nr, pt, pt_init, q, qc, qr, zu160 ONLY: hyp, pt_init, zu 154 161 155 162 USE cloud_parameters, & 156 ONLY: cp, hyrho, nc_const, pt_d_t, r_d, t_d_pt163 ONLY: cp, hyrho, prr, pt_d_t, r_d, t_d_pt 157 164 158 165 USE control_parameters, & 159 166 ONLY: call_microphysics_at_all_substeps, drizzle, dt_3d, dt_micro, & 160 g, intermediate_timestep_count, 161 l arge_scale_forcing, lsf_surf, precipitation, pt_surface,&162 rho_surface,surface_pressure167 g, intermediate_timestep_count, large_scale_forcing, & 168 lsf_surf, microphysics_kessler, microphysics_seifert, & 169 pt_surface, rho_surface,surface_pressure 163 170 164 171 USE indices, & … … 172 179 IMPLICIT NONE 173 180 174 INTEGER(iwp) :: i !<175 INTEGER(iwp) :: j !<176 181 INTEGER(iwp) :: k !< 177 182 … … 193 198 hyrho(k) = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) ) 194 199 ENDDO 200 195 201 ! 196 202 !-- Compute reference density … … 207 213 208 214 ! 215 !-- Reset precipitation rate 216 IF ( intermediate_timestep_count == 1 ) prr = 0.0_wp 217 218 ! 209 219 !-- Compute cloud physics 210 IF ( precipitation ) THEN 220 IF ( microphysics_kessler ) THEN 221 222 CALL autoconversion_kessler 223 224 ELSEIF ( microphysics_seifert ) THEN 225 211 226 CALL adjust_cloud 212 227 CALL autoconversion … … 215 230 CALL evaporation_rain 216 231 CALL sedimentation_rain 232 233 IF ( drizzle ) CALL sedimentation_cloud 234 217 235 ENDIF 218 236 219 IF ( drizzle ) CALL sedimentation_cloud 220 221 IF ( precipitation ) THEN 222 CALL calc_precipitation_amount 223 ENDIF 237 CALL calc_precipitation_amount 224 238 225 239 END SUBROUTINE microphysics_control … … 238 252 239 253 USE cloud_parameters, & 240 ONLY: eps_sb, xrmin, xrmax, hyrho , k_cc, x0254 ONLY: eps_sb, xrmin, xrmax, hyrho 241 255 242 256 USE cpulog, & … … 244 258 245 259 USE indices, & 246 ONLY: nxl, nxr, nys, nyn, nzb , nzb_s_inner, nzt260 ONLY: nxl, nxr, nys, nyn, nzb_s_inner, nzt 247 261 248 262 USE kinds … … 303 317 304 318 USE indices, & 305 ONLY: nxl, nxr, nys, nyn, nzb , nzb_s_inner, nzt319 ONLY: nxl, nxr, nys, nyn, nzb_s_inner, nzt 306 320 307 321 USE kinds … … 323 337 REAL(wp) :: rc !< 324 338 REAL(wp) :: re_lambda !< 325 REAL(wp) :: selfcoll !<326 339 REAL(wp) :: sigma_cc !< 327 340 REAL(wp) :: tau_cloud !< … … 417 430 ! Description: 418 431 ! ------------ 432 !> Autoconversion process (Kessler, 1969). 433 !------------------------------------------------------------------------------! 434 SUBROUTINE autoconversion_kessler 435 436 USE arrays_3d, & 437 ONLY: dzw, pt, q, qc 438 439 USE cloud_parameters, & 440 ONLY: l_d_cp, pt_d_t, prec_time_const, prr, ql_crit 441 442 USE control_parameters, & 443 ONLY: dt_micro 444 445 USE indices, & 446 ONLY: nxl, nxr, nyn, nys, nzb_2d, nzt 447 448 USE kinds 449 450 451 IMPLICIT NONE 452 453 INTEGER(iwp) :: i !< 454 INTEGER(iwp) :: j !< 455 INTEGER(iwp) :: k !< 456 457 REAL(wp) :: dqdt_precip !< 458 459 DO i = nxl, nxr 460 DO j = nys, nyn 461 DO k = nzb_2d(j,i)+1, nzt 462 463 IF ( qc(k,j,i) > ql_crit ) THEN 464 dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit ) 465 ELSE 466 dqdt_precip = 0.0_wp 467 ENDIF 468 469 qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro 470 q(k,j,i) = q(k,j,i) - dqdt_precip * dt_micro 471 pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * l_d_cp * pt_d_t(k) 472 473 ! 474 !-- Compute the rain rate 475 prr(nzb_2d(j,i)+1,j,i) = prr(nzb_2d(j,i)+1,j,i) + dqdt_precip * dzw(k) 476 477 ENDDO 478 ENDDO 479 ENDDO 480 481 END SUBROUTINE autoconversion_kessler 482 483 484 !------------------------------------------------------------------------------! 485 ! Description: 486 ! ------------ 419 487 !> Accretion rate (Seifert and Beheng, 2006). 420 488 !------------------------------------------------------------------------------! … … 434 502 435 503 USE indices, & 436 ONLY: nxl, nxr, nys, nyn, nzb , nzb_s_inner, nzt504 ONLY: nxl, nxr, nys, nyn, nzb_s_inner, nzt 437 505 438 506 USE kinds … … 448 516 REAL(wp) :: phi_ac !< 449 517 REAL(wp) :: tau_cloud !< 450 REAL(wp) :: xc !<451 518 452 519 CALL cpu_log( log_point_s(56), 'accretion', 'start' ) … … 516 583 517 584 USE indices, & 518 ONLY: nxl, nxr, nys, nyn, nzb , nzb_s_inner, nzt585 ONLY: nxl, nxr, nys, nyn, nzb_s_inner, nzt 519 586 520 587 USE kinds … … 579 646 580 647 USE cloud_parameters, & 581 ONLY: a_term, a_vent, b_term, b_vent, c_evap, c_term, diff_coeff_l,&582 d pirho_l, eps_sb, hyrho, kin_vis_air, k_st, l_d_cp, l_d_r,&583 l_ v, rho_l, r_v, schmidt_p_1d3, thermal_conductivity_l,&584 t _d_pt, ventilation_effect648 ONLY: a_term, a_vent, b_term, b_vent, c_evap, c_term, & 649 diff_coeff_l, dpirho_l, eps_sb, hyrho, kin_vis_air, & 650 l_d_cp, l_d_r, l_v, r_v, schmidt_p_1d3, & 651 thermal_conductivity_l, t_d_pt, ventilation_effect 585 652 586 653 USE constants, & … … 594 661 595 662 USE indices, & 596 ONLY: nxl, nxr, nys, nyn, nzb , nzb_s_inner, nzt663 ONLY: nxl, nxr, nys, nyn, nzb_s_inner, nzt 597 664 598 665 USE kinds … … 742 809 ONLY: eps_sb, hyrho, l_d_cp, nc_const, prr, pt_d_t, sed_qc_const 743 810 744 USE constants, &745 ONLY: pi746 747 811 USE control_parameters, & 748 812 ONLY: call_microphysics_at_all_substeps, dt_micro, & 749 intermediate_timestep_count , precipitation813 intermediate_timestep_count 750 814 751 815 USE cpulog, & … … 798 862 ! 799 863 !-- Compute the precipitation rate due to cloud (fog) droplets 800 IF ( precipitation ) THEN 801 IF ( call_microphysics_at_all_substeps ) THEN 802 prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) & 803 * weight_substep(intermediate_timestep_count) 804 ELSE 805 prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) 806 ENDIF 864 IF ( call_microphysics_at_all_substeps ) THEN 865 prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) & 866 * weight_substep(intermediate_timestep_count) 867 ELSE 868 prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) 807 869 ENDIF 808 870 … … 828 890 829 891 USE cloud_parameters, & 830 ONLY: a_term, b_term, c_term, cof, dpirho_l, eps_sb, hyrho,&831 limiter_sedimentation, l_d_cp, prr, pt_d_t, stp892 ONLY: a_term, b_term, c_term, dpirho_l, eps_sb, hyrho, & 893 limiter_sedimentation, l_d_cp, prr, pt_d_t 832 894 833 895 USE control_parameters, & … … 857 919 REAL(wp) :: d_min !< 858 920 REAL(wp) :: dr !< 859 REAL(wp) :: dt_sedi !<860 921 REAL(wp) :: flux !< 861 922 REAL(wp) :: lambda_r !< … … 865 926 REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !< 866 927 REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !< 867 REAL(wp), DIMENSION(nzb:nzt+1) :: d_nr !<868 REAL(wp), DIMENSION(nzb:nzt+1) :: d_qr !<869 928 REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !< 870 929 REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !< … … 876 935 CALL cpu_log( log_point_s(60), 'sed_rain', 'start' ) 877 936 878 IF ( intermediate_timestep_count == 1 ) prr(:,:,:) = 0.0_wp879 937 ! 880 938 !-- Compute velocities … … 1105 1163 1106 1164 USE cloud_parameters, & 1107 ONLY: cp, hyrho, nc_const, p t_d_t, r_d, t_d_pt1165 ONLY: cp, hyrho, nc_const, prr, pt_d_t, r_d, t_d_pt 1108 1166 1109 1167 USE control_parameters, & 1110 1168 ONLY: call_microphysics_at_all_substeps, drizzle, dt_3d, dt_micro, & 1111 1169 g, intermediate_timestep_count, large_scale_forcing, & 1112 lsf_surf, precipitation, pt_surface,&1113 rho_surface,surface_pressure1170 lsf_surf, microphysics_seifert, microphysics_kessler, & 1171 pt_surface, rho_surface, surface_pressure 1114 1172 1115 1173 USE indices, & … … 1162 1220 qc_1d(:) = qc(:,j,i) 1163 1221 nc_1d(:) = nc_const 1164 IF ( precipitation) THEN1222 IF ( microphysics_seifert ) THEN 1165 1223 qr_1d(:) = qr(:,j,i) 1166 1224 nr_1d(:) = nr(:,j,i) … … 1168 1226 1169 1227 ! 1228 !-- Reset precipitation rate 1229 IF ( intermediate_timestep_count == 1 ) prr(:,j,i) = 0.0_wp 1230 1231 ! 1170 1232 !-- Compute cloud physics 1171 IF ( precipitation ) THEN 1172 CALL adjust_cloud( i,j ) 1233 IF( microphysics_kessler ) THEN 1234 1235 CALL autoconversion_kessler( i,j ) 1236 1237 ELSEIF ( microphysics_seifert ) THEN 1238 1239 CALL adjust_cloud( i,j ) 1173 1240 CALL autoconversion( i,j ) 1174 1241 CALL accretion( i,j ) … … 1176 1243 CALL evaporation_rain( i,j ) 1177 1244 CALL sedimentation_rain( i,j ) 1245 1246 IF ( drizzle ) CALL sedimentation_cloud( i,j ) 1247 1178 1248 ENDIF 1179 1249 1180 IF ( drizzle ) CALL sedimentation_cloud( i,j ) 1181 1182 IF ( precipitation ) THEN 1183 CALL calc_precipitation_amount( i,j ) 1184 ENDIF 1250 CALL calc_precipitation_amount( i,j ) 1185 1251 1186 1252 ! … … 1188 1254 q(:,j,i) = q_1d(:) 1189 1255 pt(:,j,i) = pt_1d(:) 1190 IF ( precipitation) THEN1256 IF ( microphysics_seifert ) THEN 1191 1257 qr(:,j,i) = qr_1d(:) 1192 1258 nr(:,j,i) = nr_1d(:) … … 1210 1276 1211 1277 USE cloud_parameters, & 1212 ONLY: eps_sb, xrmin, xrmax, hyrho , k_cc, x01278 ONLY: eps_sb, xrmin, xrmax, hyrho 1213 1279 1214 1280 USE indices, & 1215 ONLY: nzb , nzb_s_inner, nzt1281 ONLY: nzb_s_inner, nzt 1216 1282 1217 1283 USE kinds … … 1258 1324 USE cloud_parameters, & 1259 1325 ONLY: a_1, a_2, a_3, b_1, b_2, b_3, beta_cc, c_1, c_2, c_3, & 1260 c_const, dpirho_l, eps_sb, hyrho, k _cc, kin_vis_air, x01326 c_const, dpirho_l, eps_sb, hyrho, kin_vis_air, k_cc, x0 1261 1327 1262 1328 USE control_parameters, & … … 1267 1333 1268 1334 USE indices, & 1269 ONLY: nzb , nzb_s_inner, nzt1335 ONLY: nzb_s_inner, nzt 1270 1336 1271 1337 USE kinds … … 1287 1353 REAL(wp) :: rc !< 1288 1354 REAL(wp) :: re_lambda !< 1289 REAL(wp) :: selfcoll !<1290 1355 REAL(wp) :: sigma_cc !< 1291 1356 REAL(wp) :: tau_cloud !< … … 1366 1431 END SUBROUTINE autoconversion_ij 1367 1432 1433 !------------------------------------------------------------------------------! 1434 ! Description: 1435 ! ------------ 1436 !> Autoconversion process (Kessler, 1969). 1437 !------------------------------------------------------------------------------! 1438 SUBROUTINE autoconversion_kessler_ij( i, j ) 1439 1440 USE arrays_3d, & 1441 ONLY: dzw, pt_1d, q_1d, qc_1d 1442 1443 USE cloud_parameters, & 1444 ONLY: l_d_cp, pt_d_t, prec_time_const, prr, ql_crit 1445 1446 USE control_parameters, & 1447 ONLY: dt_micro 1448 1449 USE indices, & 1450 ONLY: nzb_2d, nzt 1451 1452 USE kinds 1453 1454 1455 IMPLICIT NONE 1456 1457 INTEGER(iwp) :: i !< 1458 INTEGER(iwp) :: j !< 1459 INTEGER(iwp) :: k !< 1460 1461 REAL(wp) :: dqdt_precip !< 1462 1463 DO k = nzb_2d(j,i)+1, nzt 1464 1465 IF ( qc_1d(k) > ql_crit ) THEN 1466 dqdt_precip = prec_time_const * ( qc_1d(k) - ql_crit ) 1467 ELSE 1468 dqdt_precip = 0.0_wp 1469 ENDIF 1470 1471 qc_1d(k) = qc_1d(k) - dqdt_precip * dt_micro 1472 q_1d(k) = q_1d(k) - dqdt_precip * dt_micro 1473 pt_1d(k) = pt_1d(k) + dqdt_precip * dt_micro * l_d_cp * pt_d_t(k) 1474 1475 ! 1476 !-- Compute the rain rate 1477 prr(nzb_2d(j,i)+1,j,i) = prr(nzb_2d(j,i)+1,j,i) + & 1478 dqdt_precip * dzw(k) 1479 1480 ENDDO 1481 1482 END SUBROUTINE autoconversion_kessler_ij 1368 1483 1369 1484 !------------------------------------------------------------------------------! … … 1384 1499 1385 1500 USE indices, & 1386 ONLY: nzb , nzb_s_inner, nzt1501 ONLY: nzb_s_inner, nzt 1387 1502 1388 1503 USE kinds … … 1398 1513 REAL(wp) :: phi_ac !< 1399 1514 REAL(wp) :: tau_cloud !< 1400 REAL(wp) :: xc !<1401 1515 1402 1516 DO k = nzb_s_inner(j,i)+1, nzt … … 1453 1567 1454 1568 USE indices, & 1455 ONLY: nzb , nzb_s_inner, nzt1569 ONLY: nzb_s_inner, nzt 1456 1570 1457 1571 USE kinds … … 1507 1621 USE cloud_parameters, & 1508 1622 ONLY: a_term, a_vent, b_term, b_vent, c_evap, c_term, diff_coeff_l,& 1509 dpirho_l, eps_sb, hyrho, kin_vis_air, k_st, l_d_cp, l_d_r,&1510 l_v, r ho_l, r_v, schmidt_p_1d3, thermal_conductivity_l,&1623 dpirho_l, eps_sb, hyrho, kin_vis_air, l_d_cp, l_d_r, & 1624 l_v, r_v, schmidt_p_1d3, thermal_conductivity_l, & 1511 1625 t_d_pt, ventilation_effect 1512 1626 … … 1518 1632 1519 1633 USE indices, & 1520 ONLY: nzb , nzb_s_inner, nzt1634 ONLY: nzb_s_inner, nzt 1521 1635 1522 1636 USE kinds … … 1655 1769 ONLY: eps_sb, hyrho, l_d_cp, prr, pt_d_t, sed_qc_const 1656 1770 1657 USE constants, &1658 ONLY: pi1659 1660 1771 USE control_parameters, & 1661 ONLY: call_microphysics_at_all_substeps, dt_ do2d_xy, dt_micro,&1662 intermediate_timestep_count , precipitation1772 ONLY: call_microphysics_at_all_substeps, dt_micro, & 1773 intermediate_timestep_count 1663 1774 1664 1775 USE indices, & … … 1701 1812 ! 1702 1813 !-- Compute the precipitation rate of cloud (fog) droplets 1703 IF ( precipitation ) THEN 1704 IF ( call_microphysics_at_all_substeps ) THEN 1705 prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * & 1814 IF ( call_microphysics_at_all_substeps ) THEN 1815 prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * & 1706 1816 weight_substep(intermediate_timestep_count) 1707 ELSE 1708 prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) 1709 ENDIF 1817 ELSE 1818 prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) 1710 1819 ENDIF 1711 1820 … … 1727 1836 1728 1837 USE cloud_parameters, & 1729 ONLY: a_term, b_term, c_term, cof, dpirho_l, eps_sb, hyrho, & 1730 limiter_sedimentation, l_d_cp, precipitation_amount, prr, & 1731 pt_d_t, stp 1838 ONLY: a_term, b_term, c_term, dpirho_l, eps_sb, hyrho, & 1839 limiter_sedimentation, l_d_cp, prr, pt_d_t 1732 1840 1733 1841 USE control_parameters, & 1734 ONLY: call_microphysics_at_all_substeps, dt_do2d_xy, dt_micro, & 1735 dt_3d, intermediate_timestep_count, & 1736 intermediate_timestep_count_max, & 1737 precipitation_amount_interval, time_do2d_xy 1842 ONLY: call_microphysics_at_all_substeps, dt_micro, & 1843 intermediate_timestep_count 1738 1844 1739 1845 USE indices, & … … 1757 1863 REAL(wp) :: d_min !< 1758 1864 REAL(wp) :: dr !< 1759 REAL(wp) :: dt_sedi !<1760 1865 REAL(wp) :: flux !< 1761 1866 REAL(wp) :: lambda_r !< … … 1765 1870 REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !< 1766 1871 REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !< 1767 REAL(wp), DIMENSION(nzb:nzt+1) :: d_nr !<1768 REAL(wp), DIMENSION(nzb:nzt+1) :: d_qr !<1769 1872 REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !< 1770 1873 REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !< … … 1774 1877 REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr !< 1775 1878 1776 IF ( intermediate_timestep_count == 1 ) prr(:,j,i) = 0.0_wp1777 1879 ! 1778 1880 !-- Compute velocities … … 1939 2041 ONLY: call_microphysics_at_all_substeps, dt_do2d_xy, dt_3d, & 1940 2042 intermediate_timestep_count, intermediate_timestep_count_max,& 1941 precipitation_amount_interval, & 1942 time_do2d_xy 2043 precipitation_amount_interval, time_do2d_xy 1943 2044 1944 2045 USE indices, &
Note: See TracChangeset
for help on using the changeset viewer.
