Changeset 2292 for palm/trunk/SOURCE/microphysics_mod.f90
- Timestamp:
- Jun 20, 2017 9:51:42 AM (7 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/microphysics_mod.f90 (modified) (51 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/microphysics_mod.f90
r2233 r2292 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 28 ! includes two more prognostic equations for cloud drop concentration (nc) 29 ! and cloud water content (qc). 30 ! - The process of activation is parameterized with a simple Twomey 31 ! activion scheme or with considering solution and curvature 32 ! effects (Khvorostyanov and Curry ,2006). 33 ! - The saturation adjustment scheme is replaced by the parameterization 34 ! of condensation rates (Khairoutdinov and Kogan, 2000, Mon. Wea. Rev.,128). 35 ! - All other microphysical processes of Seifert and Beheng are used. 36 ! Additionally, in those processes the reduction of cloud number concentration 37 ! is considered. 38 ! 39 ! 2233 2017-05-30 18:08:54Z suehring 27 40 ! 28 41 ! 2232 2017-05-30 17:47:52Z suehring … … 123 136 IMPLICIT NONE 124 137 125 LOGICAL :: cloud_water_sedimentation = .FALSE. !< cloud water sedimentation 126 LOGICAL :: limiter_sedimentation = .TRUE. !< sedimentation limiter 127 LOGICAL :: collision_turbulence = .FALSE. !< turbulence effects 128 LOGICAL :: ventilation_effect = .TRUE. !< ventilation effect 138 LOGICAL :: cloud_water_sedimentation = .FALSE. !< cloud water sedimentation 139 LOGICAL :: curvature_solution_effects_bulk = .FALSE. !< flag for considering koehler theory 140 LOGICAL :: limiter_sedimentation = .TRUE. !< sedimentation limiter 141 LOGICAL :: collision_turbulence = .FALSE. !< turbulence effects 142 LOGICAL :: ventilation_effect = .TRUE. !< ventilation effect 129 143 130 144 REAL(wp) :: a_1 = 8.69E-4_wp !< coef. in turb. parametrization (cm-2 s3) … … 147 161 REAL(wp) :: diff_coeff_l = 0.23E-4_wp !< diffusivity of water vapor (m2 s-1) 148 162 REAL(wp) :: eps_sb = 1.0E-10_wp !< threshold in two-moments scheme 163 REAL(wp) :: eps_mr = 0.0_wp !< threshold for morrison scheme 149 164 REAL(wp) :: k_cc = 9.44E09_wp !< const. cloud-cloud kernel (m3 kg-2 s-1) 150 165 REAL(wp) :: k_cr0 = 4.33_wp !< const. cloud-rain kernel (m3 kg-1 s-1) … … 161 176 REAL(wp) :: w_precipitation = 9.65_wp !< maximum terminal velocity (m s-1) 162 177 REAL(wp) :: x0 = 2.6E-10_wp !< separating drop mass (kg) 178 REAL(wp) :: xamin = 5.24E-19_wp !< average aerosol mass (kg) (~ 0.05µm) 179 REAL(wp) :: xcmin = 4.18E-15_wp !< minimum cloud drop size (kg) (~ 1µm) 180 REAL(wp) :: xcmax = 2.6E-10_wp !< maximum cloud drop size (kg) (~ 40µm) 163 181 REAL(wp) :: xrmin = 2.6E-10_wp !< minimum rain drop size (kg) 164 182 REAL(wp) :: xrmax = 5.0E-6_wp !< maximum rain drop site (kg) … … 167 185 REAL(wp) :: dpirho_l !< 6.0 / ( pi * rho_l ) 168 186 REAL(wp) :: dt_micro !< microphysics time step 187 REAL(wp) :: na_init = 100.0E6_wp !< Total particle/aerosol concentration (cm-3) 169 188 REAL(wp) :: nc_const = 70.0E6_wp !< cloud droplet concentration 170 189 REAL(wp) :: dt_precipitation = 100.0_wp !< timestep precipitation (s) … … 184 203 PUBLIC microphysics_control, microphysics_init 185 204 186 PUBLIC cloud_water_sedimentation, collision_turbulence, c_sedimentation, & 187 dt_precipitation, limiter_sedimentation, nc_const, sigma_gc, & 205 PUBLIC cloud_water_sedimentation, collision_turbulence, & 206 curvature_solution_effects_bulk, c_sedimentation, dt_precipitation, & 207 limiter_sedimentation, na_init, nc_const, sigma_gc, & 188 208 ventilation_effect 209 189 210 190 211 INTERFACE microphysics_control … … 198 219 END INTERFACE adjust_cloud 199 220 221 INTERFACE activation 222 MODULE PROCEDURE activation 223 MODULE PROCEDURE activation_ij 224 END INTERFACE activation 225 226 INTERFACE condensation 227 MODULE PROCEDURE condensation 228 MODULE PROCEDURE condensation_ij 229 END INTERFACE condensation 230 200 231 INTERFACE autoconversion 201 232 MODULE PROCEDURE autoconversion … … 256 287 257 288 USE control_parameters, & 258 ONLY: microphysics_ seifert289 ONLY: microphysics_morrison, microphysics_seifert 259 290 260 291 USE indices, & … … 283 314 ! 284 315 !-- Allocate 1D microphysics arrays 285 ALLOCATE ( nc_1d(nzb:nzt+1),pt_1d(nzb:nzt+1), q_1d(nzb:nzt+1), &316 ALLOCATE ( pt_1d(nzb:nzt+1), q_1d(nzb:nzt+1), & 286 317 qc_1d(nzb:nzt+1) ) 318 319 IF ( microphysics_morrison .OR. microphysics_seifert ) THEN 320 ALLOCATE ( nc_1d(nzb:nzt+1) ) 321 ENDIF 287 322 288 323 IF ( microphysics_seifert ) THEN … … 290 325 ENDIF 291 326 292 !293 !-- Initialze nc_1d with nc_const294 nc_1d = nc_const295 296 327 END SUBROUTINE microphysics_init 297 328 … … 313 344 ONLY: call_microphysics_at_all_substeps, dt_3d, g, & 314 345 intermediate_timestep_count, large_scale_forcing, & 315 lsf_surf, microphysics_kessler, microphysics_seifert, & 316 pt_surface, rho_surface,surface_pressure 346 lsf_surf, microphysics_kessler, microphysics_morrison, & 347 microphysics_seifert, pt_surface, & 348 rho_surface,surface_pressure 317 349 318 350 USE indices, & … … 373 405 374 406 CALL adjust_cloud 407 IF ( microphysics_morrison ) CALL activation 408 IF ( microphysics_morrison ) CALL condensation 375 409 CALL autoconversion 376 410 CALL accretion … … 396 430 397 431 USE arrays_3d, & 398 ONLY: q r, nr432 ONLY: qc, nc, qr, nr 399 433 400 434 USE cloud_parameters, & 401 435 ONLY: hyrho 436 437 USE control_parameters, & 438 ONLY: microphysics_morrison 402 439 403 440 USE cpulog, & … … 434 471 ENDIF 435 472 ENDIF 473 IF ( microphysics_morrison ) THEN 474 IF ( qc(k,j,i) <= eps_sb ) THEN 475 qc(k,j,i) = 0.0_wp 476 nc(k,j,i) = 0.0_wp 477 ELSE 478 IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) ) THEN 479 nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin * & 480 MERGE( 1.0_wp, 0.0_wp, & 481 BTEST( wall_flags_0(k,j,i), 0 ) ) 482 ENDIF 483 ENDIF 484 ENDIF 436 485 ENDDO 437 486 ENDDO … … 442 491 END SUBROUTINE adjust_cloud 443 492 493 !------------------------------------------------------------------------------! 494 ! Description: 495 ! ------------ 496 !> Calculate number of activated condensation nucleii after simple activation 497 !> scheme of Twomey, 1959. 498 !------------------------------------------------------------------------------! 499 SUBROUTINE activation 500 501 USE arrays_3d, & 502 ONLY: hyp, nc, nr, pt, q, qc, qr 503 504 USE cloud_parameters, & 505 ONLY: hyrho, l_d_cp, l_d_r, l_v, rho_l, r_v, t_d_pt 506 507 USE constants, & 508 ONLY: pi 509 510 USE cpulog, & 511 ONLY: cpu_log, log_point_s 512 513 USE indices, & 514 ONLY: nxlg, nxrg, nysg, nyng, nzb, nzt 515 516 USE kinds 517 518 USE control_parameters, & 519 ONLY: simulated_time 520 521 USE particle_attributes, & 522 ONLY: molecular_weight_of_solute, molecular_weight_of_water, rho_s, & 523 s1, s2, s3, vanthoff 524 525 IMPLICIT NONE 526 527 INTEGER(iwp) :: i !< 528 INTEGER(iwp) :: j !< 529 INTEGER(iwp) :: k !< 530 531 REAL(wp) :: activ !< 532 REAL(wp) :: afactor !< 533 REAL(wp) :: alpha !< 534 REAL(wp) :: beta_act !< 535 REAL(wp) :: bfactor !< 536 REAL(wp) :: e_s !< 537 REAL(wp) :: k_act !< 538 REAL(wp) :: n_act !< 539 REAL(wp) :: n_ccn !< 540 REAL(wp) :: q_s !< 541 REAL(wp) :: rd0 !< 542 REAL(wp) :: s_0 !< 543 REAL(wp) :: sat !< 544 REAL(wp) :: sat_max !< 545 REAL(wp) :: sigma !< 546 REAL(wp) :: t_int !< 547 REAL(wp) :: t_l !< 548 549 CALL cpu_log( log_point_s(65), 'activation', 'start' ) 550 551 DO i = nxlg, nxrg 552 DO j = nysg, nyng 553 DO k = nzb+1, nzt 554 555 ! 556 !-- Actual liquid water temperature: 557 t_l = t_d_pt(k) * pt(k,j,i) 558 559 ! 560 !-- Calculate actual temperature 561 t_int = pt(k,j,i) * ( hyp(k) / 100000.0_wp )**0.286_wp 562 ! 563 !-- Saturation vapor pressure at t_l: 564 e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / & 565 ( t_l - 35.86_wp ) & 566 ) 567 ! 568 !-- Computation of saturation humidity: 569 q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s ) 570 alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l ) 571 q_s = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / & 572 ( 1.0_wp + alpha * q_s ) 573 574 !-- Supersaturation: 575 sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp 576 577 ! 578 !-- Prescribe parameters for activation 579 !-- (see: Bott + Trautmann, 2002, Atm. Res., 64) 580 k_act = 0.7_wp 581 582 IF ( sat >= 0.0 .AND. .NOT. curvature_solution_effects_bulk ) THEN 583 ! 584 !-- Compute the number of activated Aerosols 585 !-- (see: Twomey, 1959, Pure and applied Geophysics, 43) 586 n_act = na_init * sat**k_act 587 ! 588 !-- Compute the number of cloud droplets 589 !-- (see: Morrison + Grabowski, 2007, JAS, 64) 590 ! activ = MAX( n_act - nc(k,j,i), 0.0_wp) / dt_micro 591 592 ! 593 !-- Compute activation rate after Khairoutdinov and Kogan 594 !-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128) 595 sat_max = 1.0_wp / 100.0_wp 596 activ = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN & 597 ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) / & 598 dt_micro 599 ELSEIF ( sat >= 0.0 .AND. curvature_solution_effects_bulk ) THEN 600 ! 601 !-- Curvature effect (afactor) with surface tension 602 !-- parameterization by Straka (2009) 603 sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp ) 604 afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int ) 605 ! 606 !-- Solute effect (bfactor) 607 bfactor = vanthoff * molecular_weight_of_water * & 608 rho_s / ( molecular_weight_of_solute * rho_l ) 609 610 ! 611 !-- Prescribe power index that describes the soluble fraction 612 !-- of an aerosol particle (beta) and mean geometric radius of 613 !-- dry aerosol spectrum 614 !-- (see: Morrison + Grabowski, 2007, JAS, 64) 615 beta_act = 0.5_wp 616 rd0 = 0.05E-6_wp 617 ! 618 !-- Calculate mean geometric supersaturation (s_0) with 619 !-- parameterization by Khvorostyanov and Curry (2006) 620 s_0 = rd0 **(- ( 1.0_wp + beta_act ) ) * & 621 ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp 622 623 ! 624 !-- Calculate number of activated CCN as a function of 625 !-- supersaturation and taking Koehler theory into account 626 !-- (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111) 627 n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF( & 628 LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(s1) ) ) ) 629 activ = MAX( ( n_ccn - nc(k,j,i) ) / dt_micro, 0.0_wp ) 630 ENDIF 631 632 nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro), na_init) 633 634 ENDDO 635 ENDDO 636 ENDDO 637 638 CALL cpu_log( log_point_s(65), 'activation', 'stop' ) 639 640 END SUBROUTINE activation 641 642 643 !------------------------------------------------------------------------------! 644 ! Description: 645 ! ------------ 646 !> Calculate condensation rate for cloud water content (after Khairoutdinov and 647 !> Kogan, 2000). 648 !------------------------------------------------------------------------------! 649 SUBROUTINE condensation 650 651 USE arrays_3d, & 652 ONLY: hyp, nr, pt, q, qc, qr, nc 653 654 USE cloud_parameters, & 655 ONLY: hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt 656 657 USE constants, & 658 ONLY: pi 659 660 USE cpulog, & 661 ONLY: cpu_log, log_point_s 662 663 USE indices, & 664 ONLY: nxlg, nxrg, nysg, nyng, nzb, nzt 665 666 USE kinds 667 668 USE control_parameters, & 669 ONLY: simulated_time 670 671 IMPLICIT NONE 672 673 INTEGER(iwp) :: i !< 674 INTEGER(iwp) :: j !< 675 INTEGER(iwp) :: k !< 676 677 REAL(wp) :: alpha !< 678 REAL(wp) :: cond !< 679 REAL(wp) :: cond_max !< 680 REAL(wp) :: dc !< 681 REAL(wp) :: e_s !< 682 REAL(wp) :: evap !< 683 REAL(wp) :: evap_nc !< 684 REAL(wp) :: g_fac !< 685 REAL(wp) :: nc_0 !< 686 REAL(wp) :: q_s !< 687 REAL(wp) :: sat !< 688 REAL(wp) :: t_l !< 689 REAL(wp) :: temp !< 690 REAL(wp) :: xc !< 691 692 CALL cpu_log( log_point_s(66), 'condensation', 'start' ) 693 694 DO i = nxlg, nxrg 695 DO j = nysg, nyng 696 DO k = nzb+1, nzt 697 ! 698 !-- Actual liquid water temperature: 699 t_l = t_d_pt(k) * pt(k,j,i) 700 ! 701 !-- Saturation vapor pressure at t_l: 702 e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / & 703 ( t_l - 35.86_wp ) & 704 ) 705 ! 706 !-- Computation of saturation humidity: 707 q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s ) 708 alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l ) 709 q_s = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / & 710 ( 1.0_wp + alpha * q_s ) 711 712 !-- Supersaturation: 713 sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp 714 715 ! 716 !-- Actual temperature: 717 temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) ) 718 719 g_fac = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * & 720 l_v / ( thermal_conductivity_l * temp ) & 721 + r_v * temp / ( diff_coeff_l * e_s ) & 722 ) 723 ! 724 !-- Mean weight of cloud drops 725 IF ( nc(k,j,i) <= 0.0_wp) CYCLE 726 xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin) 727 ! 728 !-- Weight averaged diameter of cloud drops: 729 dc = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp ) 730 ! 731 !-- Integral diameter of cloud drops 732 nc_0 = nc(k,j,i) * dc 733 ! 734 !-- Condensation needs only to be calculated in supersaturated regions 735 IF ( sat > 0.0_wp ) THEN 736 ! 737 !-- Condensation rate of cloud water content 738 !-- after KK scheme. 739 !-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128) 740 cond = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k) 741 cond_max = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i) 742 cond = MIN( cond, cond_max / dt_micro ) 743 744 qc(k,j,i) = qc(k,j,i) + cond * dt_micro 745 ELSEIF ( sat < 0.0_wp ) THEN 746 evap = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k) 747 evap = MAX( evap, -qc(k,j,i) / dt_micro ) 748 749 qc(k,j,i) = qc(k,j,i) + evap * dt_micro 750 ENDIF 751 IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) ) THEN 752 nc(k,j,i) = qc(k,j,i) * hyrho(k) / xcmin 753 ENDIF 754 ENDDO 755 ENDDO 756 ENDDO 757 758 CALL cpu_log( log_point_s(66), 'condensation', 'stop' ) 759 760 END SUBROUTINE condensation 761 444 762 445 763 !------------------------------------------------------------------------------! … … 451 769 452 770 USE arrays_3d, & 453 ONLY: diss, dzu, n r, qc, qr771 ONLY: diss, dzu, nc, nr, qc, qr 454 772 455 773 USE cloud_parameters, & … … 457 775 458 776 USE control_parameters, & 459 ONLY: rho_surface777 ONLY: microphysics_morrison, rho_surface 460 778 461 779 USE cpulog, & … … 482 800 REAL(wp) :: k_au !< 483 801 REAL(wp) :: l_mix !< 802 REAL(wp) :: nc_auto !< 484 803 REAL(wp) :: nu_c !< 485 804 REAL(wp) :: phi_au !< … … 500 819 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 501 820 821 nc_auto = MERGE( nc(k,j,i), nc_const, microphysics_morrison ) 822 502 823 IF ( qc(k,j,i) > eps_sb ) THEN 503 824 … … 518 839 ! 519 840 !-- Mean weight of cloud droplets: 520 xc = hyrho(k) * qc(k,j,i) / nc_ const841 xc = hyrho(k) * qc(k,j,i) / nc_auto 521 842 ! 522 843 !-- Parameterized turbulence effects on autoconversion (Seifert, … … 569 890 nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro & 570 891 * flag 892 IF ( microphysics_morrison ) THEN 893 nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp * & 894 autocon / x0 * hyrho(k) * dt_micro * flag ) 895 ENDIF 571 896 572 897 ENDIF … … 654 979 655 980 USE arrays_3d, & 656 ONLY: diss, qc, qr 981 ONLY: diss, qc, qr, nc 657 982 658 983 USE cloud_parameters, & … … 660 985 661 986 USE control_parameters, & 662 ONLY: rho_surface987 ONLY: microphysics_morrison, rho_surface 663 988 664 989 USE cpulog, & … … 679 1004 REAL(wp) :: flag !< flag to mask topography grid points 680 1005 REAL(wp) :: k_cr !< 1006 REAL(wp) :: nc_accr !< 681 1007 REAL(wp) :: phi_ac !< 682 1008 REAL(wp) :: tau_cloud !< 1009 REAL(wp) :: xc !< 1010 683 1011 684 1012 CALL cpu_log( log_point_s(56), 'accretion', 'start' ) … … 691 1019 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 692 1020 693 IF ( ( qc(k,j,i) > eps_sb ) .AND. ( qr(k,j,i) > eps_sb ) ) THEN 1021 nc_accr = MERGE( nc(k,j,i), nc_const, microphysics_morrison ) 1022 1023 IF ( ( qc(k,j,i) > eps_sb ) .AND. ( qr(k,j,i) > eps_sb ) & 1024 .AND. ( nc_accr > eps_mr ) ) THEN 694 1025 ! 695 1026 !-- Intern time scale of coagulation (Seifert and Beheng, 2006): … … 699 1030 !-- Beheng, 2001): 700 1031 phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4 1032 1033 ! 1034 !-- Mean weight of cloud drops 1035 xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin) 701 1036 ! 702 1037 !-- Parameterized turbulence effects on autoconversion (Seifert, … … 719 1054 qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag 720 1055 qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag 1056 IF ( microphysics_morrison ) THEN 1057 nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), & 1058 accr / xc * hyrho(k) * dt_micro * flag) 1059 ENDIF 721 1060 722 1061 ENDIF … … 976 1315 977 1316 USE arrays_3d, & 978 ONLY: ddzu, dzu, pt, prr, q, qc1317 ONLY: ddzu, dzu, nc, pt, prr, q, qc 979 1318 980 1319 USE cloud_parameters, & … … 982 1321 983 1322 USE control_parameters, & 984 ONLY: call_microphysics_at_all_substeps, intermediate_timestep_count 1323 ONLY: call_microphysics_at_all_substeps, & 1324 intermediate_timestep_count, microphysics_morrison 985 1325 986 1326 USE cpulog, & … … 1002 1342 INTEGER(iwp) :: k !< 1003 1343 1004 REAL(wp) :: flag !< flag to mask topography grid points 1344 REAL(wp) :: flag !< flag to mask topography grid points 1345 REAL(wp) :: nc_sedi !< 1346 1005 1347 REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !< 1348 REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nc !< 1349 1006 1350 1007 1351 CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' ) … … 1015 1359 !-- Predetermine flag to mask topography 1016 1360 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 1361 nc_sedi = MERGE ( nc(k,j,i), nc_const, microphysics_morrison ) 1362 1363 ! 1364 !-- Sedimentation fluxes for number concentration are only calculated 1365 !-- for cloud_scheme = 'morrison' 1366 IF ( microphysics_morrison ) THEN 1367 IF ( qc(k,j,i) > eps_sb .AND. nc(k,j,i) > eps_mr ) THEN 1368 sed_nc(k) = sed_qc_const * & 1369 ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) * & 1370 ( nc(k,j,i) )**( 1.0_wp / 3.0_wp ) 1371 ELSE 1372 sed_nc(k) = 0.0_wp 1373 ENDIF 1374 1375 sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) * & 1376 nc(k,j,i) / dt_micro + sed_nc(k+1) & 1377 ) * flag 1378 1379 nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) * & 1380 ddzu(k+1) / hyrho(k) * dt_micro * flag 1381 ENDIF 1017 1382 1018 1383 IF ( qc(k,j,i) > eps_sb ) THEN 1019 sed_qc(k) = sed_qc_const * nc_ const**( -2.0_wp / 3.0_wp ) *&1384 sed_qc(k) = sed_qc_const * nc_sedi**( -2.0_wp / 3.0_wp ) * & 1020 1385 ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp ) * & 1021 1386 flag … … 1389 1754 1390 1755 USE arrays_3d, & 1391 ONLY: hyp, n r, pt, pt_init, prr, q, qc, qr, zu1756 ONLY: hyp, nc, nr, pt, pt_init, prr, q, qc, qr, zu 1392 1757 1393 1758 USE cloud_parameters, & … … 1397 1762 ONLY: call_microphysics_at_all_substeps, dt_3d, g, & 1398 1763 intermediate_timestep_count, large_scale_forcing, & 1399 lsf_surf, microphysics_seifert, microphysics_kessler, & 1400 pt_surface, rho_surface, surface_pressure 1764 lsf_surf, microphysics_morrison, microphysics_seifert, & 1765 microphysics_kessler, pt_surface, rho_surface, & 1766 surface_pressure 1401 1767 1402 1768 USE indices, & … … 1448 1814 pt_1d(:) = pt(:,j,i) 1449 1815 qc_1d(:) = qc(:,j,i) 1450 nc_1d(:) = nc_const1451 1816 IF ( microphysics_seifert ) THEN 1452 1817 qr_1d(:) = qr(:,j,i) 1453 1818 nr_1d(:) = nr(:,j,i) 1454 1819 ENDIF 1820 IF ( microphysics_morrison ) THEN 1821 nc_1d(:) = nc(:,j,i) 1822 ENDIF 1823 1455 1824 1456 1825 ! … … 1468 1837 1469 1838 CALL adjust_cloud( i,j ) 1839 IF ( microphysics_morrison ) CALL activation( i,j ) 1840 IF ( microphysics_morrison ) CALL condensation( i,j ) 1470 1841 CALL autoconversion( i,j ) 1471 1842 CALL accretion( i,j ) … … 1483 1854 q(:,j,i) = q_1d(:) 1484 1855 pt(:,j,i) = pt_1d(:) 1856 IF ( microphysics_morrison ) THEN 1857 qc(:,j,i) = qc_1d(:) 1858 nc(:,j,i) = nc_1d(:) 1859 ENDIF 1485 1860 IF ( microphysics_seifert ) THEN 1486 1861 qr(:,j,i) = qr_1d(:) … … 1503 1878 USE cloud_parameters, & 1504 1879 ONLY: hyrho 1880 1881 USE control_parameters, & 1882 ONLY: microphysics_morrison 1505 1883 1506 1884 USE indices, & … … 1538 1916 ENDIF 1539 1917 1918 IF ( microphysics_morrison ) THEN 1919 IF ( qc_1d(k) <= eps_sb ) THEN 1920 qc_1d(k) = 0.0_wp 1921 nc_1d(k) = 0.0_wp 1922 ELSE 1923 IF ( nc_1d(k) * xcmin > qc_1d(k) * hyrho(k) ) THEN 1924 nc_1d(k) = qc_1d(k) * hyrho(k) / xamin * flag 1925 ENDIF 1926 ENDIF 1927 ENDIF 1928 1540 1929 ENDDO 1541 1930 1542 1931 END SUBROUTINE adjust_cloud_ij 1932 1933 !------------------------------------------------------------------------------! 1934 ! Description: 1935 ! ------------ 1936 !> Calculate number of activated condensation nucleii after simple activation 1937 !> scheme of Twomey, 1959. 1938 !------------------------------------------------------------------------------! 1939 SUBROUTINE activation_ij( i, j ) 1940 1941 USE arrays_3d, & 1942 ONLY: hyp, nr, pt, q, qc, qr, nc 1943 1944 USE cloud_parameters, & 1945 ONLY: hyrho, l_d_cp, l_d_r, l_v, rho_l, r_v, t_d_pt 1946 1947 USE constants, & 1948 ONLY: pi 1949 1950 USE cpulog, & 1951 ONLY: cpu_log, log_point_s 1952 1953 USE indices, & 1954 ONLY: nxlg, nxrg, nysg, nyng, nzb, nzt 1955 1956 USE kinds 1957 1958 USE control_parameters, & 1959 ONLY: simulated_time 1960 1961 USE particle_attributes, & 1962 ONLY: molecular_weight_of_solute, molecular_weight_of_water, rho_s, & 1963 s1, s2, s3, vanthoff 1964 1965 IMPLICIT NONE 1966 1967 INTEGER(iwp) :: i !< 1968 INTEGER(iwp) :: j !< 1969 INTEGER(iwp) :: k !< 1970 1971 REAL(wp) :: activ !< 1972 REAL(wp) :: afactor !< 1973 REAL(wp) :: alpha !< 1974 REAL(wp) :: beta_act !< 1975 REAL(wp) :: bfactor !< 1976 REAL(wp) :: e_s !< 1977 REAL(wp) :: k_act !< 1978 REAL(wp) :: n_act !< 1979 REAL(wp) :: n_ccn !< 1980 REAL(wp) :: q_s !< 1981 REAL(wp) :: rd0 !< 1982 REAL(wp) :: s_0 !< 1983 REAL(wp) :: sat !< 1984 REAL(wp) :: sat_max !< 1985 REAL(wp) :: sigma !< 1986 REAL(wp) :: t_int !< 1987 REAL(wp) :: t_l !< 1988 1989 DO k = nzb+1, nzt 1990 ! 1991 !-- Actual liquid water temperature: 1992 t_l = t_d_pt(k) * pt_1d(k) 1993 1994 ! 1995 !-- Calculate actual temperature 1996 t_int = pt_1d(k) * ( hyp(k) / 100000.0_wp )**0.286_wp 1997 ! 1998 !-- Saturation vapor pressure at t_l: 1999 e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / & 2000 ( t_l - 35.86_wp ) & 2001 ) 2002 ! 2003 !-- Computation of saturation humidity: 2004 q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s ) 2005 alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l ) 2006 q_s = q_s * ( 1.0_wp + alpha * q_1d(k) ) / & 2007 ( 1.0_wp + alpha * q_s ) 2008 2009 !-- Supersaturation: 2010 sat = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp 2011 2012 ! 2013 !-- Prescribe parameters for activation 2014 !-- (see: Bott + Trautmann, 2002, Atm. Res., 64) 2015 k_act = 0.7_wp 2016 2017 IF ( sat >= 0.0 .AND. .NOT. curvature_solution_effects_bulk ) THEN 2018 ! 2019 !-- Compute the number of activated Aerosols 2020 !-- (see: Twomey, 1959, Pure and applied Geophysics, 43) 2021 n_act = na_init * sat**k_act 2022 ! 2023 !-- Compute the number of cloud droplets 2024 !-- (see: Morrison + Grabowski, 2007, JAS, 64) 2025 ! activ = MAX( n_act - nc_d1(k), 0.0_wp) / dt_micro 2026 2027 ! 2028 !-- Compute activation rate after Khairoutdinov and Kogan 2029 !-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128) 2030 sat_max = 0.8_wp / 100.0_wp 2031 activ = MAX( 0.0_wp, ( (na_init + nc_1d(k) ) * MIN & 2032 ( 1.0_wp, ( sat / sat_max )**k_act) - nc_1d(k) ) ) / & 2033 dt_micro 2034 2035 nc_1d(k) = MIN( (nc_1d(k) + activ * dt_micro), na_init) 2036 ELSEIF ( sat >= 0.0 .AND. curvature_solution_effects_bulk ) THEN 2037 ! 2038 !-- Curvature effect (afactor) with surface tension 2039 !-- parameterization by Straka (2009) 2040 sigma = 0.0761_wp - 0.000155_wp * ( t_int - 273.15_wp ) 2041 afactor = 2.0_wp * sigma / ( rho_l * r_v * t_int ) 2042 ! 2043 !-- Solute effect (bfactor) 2044 bfactor = vanthoff * molecular_weight_of_water * & 2045 rho_s / ( molecular_weight_of_solute * rho_l ) 2046 2047 ! 2048 !-- Prescribe power index that describes the soluble fraction 2049 !-- of an aerosol particle (beta) and mean geometric radius of 2050 !-- dry aerosol spectrum 2051 !-- (see: Morrison + Grabowski, 2007, JAS, 64) 2052 beta_act = 0.5_wp 2053 rd0 = 0.05E-6_wp 2054 ! 2055 !-- Calculate mean geometric supersaturation (s_0) with 2056 !-- parameterization by Khvorostyanov and Curry (2006) 2057 s_0 = rd0 **(- ( 1.0_wp + beta_act ) ) * & 2058 ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp 2059 2060 ! 2061 !-- Calculate number of activated CCN as a function of 2062 !-- supersaturation and taking Koehler theory into account 2063 !-- (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111) 2064 n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF( & 2065 LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(s1) ) ) ) 2066 activ = MAX( ( n_ccn ) / dt_micro, 0.0_wp ) 2067 2068 nc_1d(k) = MIN( (nc_1d(k) + activ * dt_micro), na_init) 2069 ENDIF 2070 2071 ENDDO 2072 2073 END SUBROUTINE activation_ij 2074 2075 !------------------------------------------------------------------------------! 2076 ! Description: 2077 ! ------------ 2078 !> Calculate condensation rate for cloud water content (after Khairoutdinov and 2079 !> Kogan, 2000). 2080 !------------------------------------------------------------------------------! 2081 SUBROUTINE condensation_ij( i, j ) 2082 2083 USE arrays_3d, & 2084 ONLY: hyp, nr, pt, q, qc, qr, nc 2085 2086 USE cloud_parameters, & 2087 ONLY: hyrho, l_d_cp, l_d_r, l_v, r_v, t_d_pt 2088 2089 USE constants, & 2090 ONLY: pi 2091 2092 USE cpulog, & 2093 ONLY: cpu_log, log_point_s 2094 2095 USE indices, & 2096 ONLY: nxlg, nxrg, nysg, nyng, nzb, nzt 2097 2098 USE kinds 2099 2100 USE control_parameters, & 2101 ONLY: simulated_time 2102 2103 IMPLICIT NONE 2104 2105 INTEGER(iwp) :: i !< 2106 INTEGER(iwp) :: j !< 2107 INTEGER(iwp) :: k !< 2108 2109 REAL(wp) :: alpha !< 2110 REAL(wp) :: cond !< 2111 REAL(wp) :: cond_max !< 2112 REAL(wp) :: dc !< 2113 REAL(wp) :: e_s !< 2114 REAL(wp) :: evap !< 2115 REAL(wp) :: evap_nc !< 2116 REAL(wp) :: g_fac !< 2117 REAL(wp) :: nc_0 !< 2118 REAL(wp) :: q_s !< 2119 REAL(wp) :: sat !< 2120 REAL(wp) :: t_l !< 2121 REAL(wp) :: temp !< 2122 REAL(wp) :: xc !< 2123 2124 2125 DO k = nzb+1, nzt 2126 ! 2127 !-- Actual liquid water temperature: 2128 t_l = t_d_pt(k) * pt_1d(k) 2129 ! 2130 !-- Saturation vapor pressure at t_l: 2131 e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / & 2132 ( t_l - 35.86_wp ) & 2133 ) 2134 ! 2135 !-- Computation of saturation humidity: 2136 q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s ) 2137 alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l ) 2138 q_s = q_s * ( 1.0_wp + alpha * q_1d(k) ) / & 2139 ( 1.0_wp + alpha * q_s ) 2140 2141 !-- Supersaturation: 2142 sat = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp 2143 2144 2145 ! 2146 !-- Actual temperature: 2147 temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) ) 2148 2149 g_fac = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * & 2150 l_v / ( thermal_conductivity_l * temp ) & 2151 + r_v * temp / ( diff_coeff_l * e_s ) & 2152 ) 2153 ! 2154 !-- Mean weight of cloud drops 2155 IF ( nc_1d(k) <= 0.0_wp) CYCLE 2156 xc = MAX( (hyrho(k) * qc_1d(k) / nc_1d(k)), xcmin) 2157 ! 2158 !-- Weight averaged diameter of cloud drops: 2159 dc = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp ) 2160 ! 2161 !-- Integral diameter of cloud drops 2162 nc_0 = nc_1d(k) * dc 2163 ! 2164 !-- Condensation needs only to be calculated in supersaturated regions 2165 IF ( sat > 0.0_wp ) THEN 2166 ! 2167 !-- Condensation rate of cloud water content 2168 !-- after KK scheme. 2169 !-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128) 2170 cond = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k) 2171 cond_max = q_1d(k) - q_s - qc_1d(k) - qr_1d(k) 2172 cond = MIN( cond, cond_max / dt_micro ) 2173 2174 qc_1d(k) = qc_1d(k) + cond * dt_micro 2175 ELSEIF ( sat < 0.0_wp ) THEN 2176 evap = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k) 2177 evap = MAX( evap, -qc_1d(k) / dt_micro ) 2178 2179 qc_1d(k) = qc_1d(k) + evap * dt_micro 2180 ENDIF 2181 ENDDO 2182 2183 END SUBROUTINE condensation_ij 1543 2184 1544 2185 … … 1557 2198 1558 2199 USE control_parameters, & 1559 ONLY: rho_surface2200 ONLY: microphysics_morrison, rho_surface 1560 2201 1561 2202 USE grid_variables, & … … 1579 2220 REAL(wp) :: k_au !< 1580 2221 REAL(wp) :: l_mix !< 2222 REAL(wp) :: nc_auto !< 1581 2223 REAL(wp) :: nu_c !< 1582 2224 REAL(wp) :: phi_au !< … … 1592 2234 !-- Predetermine flag to mask topography 1593 2235 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 2236 nc_auto = MERGE ( nc_1d(k), nc_const, microphysics_morrison ) 1594 2237 1595 2238 IF ( qc_1d(k) > eps_sb ) THEN … … 1610 2253 ! 1611 2254 !-- Mean weight of cloud droplets: 1612 xc = hyrho(k) * qc_1d(k) / nc_ 1d(k)2255 xc = hyrho(k) * qc_1d(k) / nc_auto 1613 2256 ! 1614 2257 !-- Parameterized turbulence effects on autoconversion (Seifert, … … 1658 2301 qc_1d(k) = qc_1d(k) - autocon * dt_micro * flag 1659 2302 nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro * flag 2303 IF ( microphysics_morrison ) THEN 2304 nc_1d(k) = nc_1d(k) - MIN( nc_1d(k), 2.0_wp * & 2305 autocon / x0 * hyrho(k) * dt_micro * flag ) 2306 ENDIF 1660 2307 1661 2308 ENDIF … … 1738 2385 1739 2386 USE control_parameters, & 1740 ONLY: rho_surface2387 ONLY: microphysics_morrison, rho_surface 1741 2388 1742 2389 USE indices, & … … 1754 2401 REAL(wp) :: flag !< flag to indicate first grid level above surface 1755 2402 REAL(wp) :: k_cr !< 2403 REAL(wp) :: nc_accr !< 1756 2404 REAL(wp) :: phi_ac !< 1757 2405 REAL(wp) :: tau_cloud !< 2406 REAL(wp) :: xc !< 2407 1758 2408 1759 2409 DO k = nzb+1, nzt … … 1761 2411 !-- Predetermine flag to mask topography 1762 2412 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 1763 1764 IF ( ( qc_1d(k) > eps_sb ) .AND. ( qr_1d(k) > eps_sb ) ) THEN 2413 nc_accr = MERGE ( nc_1d(k), nc_const, microphysics_morrison ) 2414 2415 IF ( ( qc_1d(k) > eps_sb ) .AND. ( qr_1d(k) > eps_sb ) .AND. & 2416 ( nc_accr > eps_mr ) ) THEN 1765 2417 ! 1766 2418 !-- Intern time scale of coagulation (Seifert and Beheng, 2006): … … 1770 2422 !-- (Seifert and Beheng, 2001): 1771 2423 phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4 2424 2425 ! 2426 !-- Mean weight of cloud drops 2427 xc = MAX( (hyrho(k) * qc_1d(k) / nc_accr), xcmin) 1772 2428 ! 1773 2429 !-- Parameterized turbulence effects on autoconversion (Seifert, … … 1784 2440 ! 1785 2441 !-- Accretion rate (Seifert and Beheng, 2006): 1786 accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac * SQRT( rho_surface * hyrho(k) ) 2442 accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac * & 2443 SQRT( rho_surface * hyrho(k) ) 1787 2444 accr = MIN( accr, qc_1d(k) / dt_micro ) 1788 2445 1789 2446 qr_1d(k) = qr_1d(k) + accr * dt_micro * flag 1790 2447 qc_1d(k) = qc_1d(k) - accr * dt_micro * flag 2448 IF ( microphysics_morrison ) THEN 2449 nc_1d(k) = nc_1d(k) - MIN( nc_1d(k), accr / xc * & 2450 hyrho(k) * dt_micro * flag & 2451 ) 2452 ENDIF 2453 1791 2454 1792 2455 ENDIF … … 2018 2681 2019 2682 USE control_parameters, & 2020 ONLY: call_microphysics_at_all_substeps, intermediate_timestep_count 2683 ONLY: call_microphysics_at_all_substeps, & 2684 intermediate_timestep_count, microphysics_morrison 2021 2685 2022 2686 USE indices, & … … 2034 2698 INTEGER(iwp) :: k !< 2035 2699 2036 REAL(wp) :: flag !< flag to indicate first grid level above surface 2700 REAL(wp) :: flag !< flag to indicate first grid level above surface 2701 REAL(wp) :: nc_sedi !< 2702 2703 REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nc !< 2037 2704 REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !< 2038 2705 … … 2043 2710 !-- Predetermine flag to mask topography 2044 2711 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 2712 nc_sedi = MERGE( nc_1d(k), nc_const, microphysics_morrison ) 2713 ! 2714 !-- Sedimentation fluxes for number concentration are only calculated 2715 !-- for cloud_scheme = 'morrison' 2716 IF ( microphysics_morrison ) THEN 2717 IF ( qc_1d(k) > eps_sb .AND. nc_1d(k) > eps_mr ) THEN 2718 sed_nc(k) = sed_qc_const * & 2719 ( qc_1d(k) * hyrho(k) )**( 2.0_wp / 3.0_wp ) * & 2720 ( nc_1d(k) )**( 1.0_wp / 3.0_wp ) 2721 ELSE 2722 sed_nc(k) = 0.0_wp 2723 ENDIF 2724 2725 sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) * & 2726 nc_1d(k) / dt_micro + sed_nc(k+1) & 2727 ) * flag 2728 2729 nc_1d(k) = nc_1d(k) + ( sed_nc(k+1) - sed_nc(k) ) * & 2730 ddzu(k+1) / hyrho(k) * dt_micro * flag 2731 ENDIF 2045 2732 2046 2733 IF ( qc_1d(k) > eps_sb ) THEN 2047 sed_qc(k) = sed_qc_const * nc_ 1d(k)**( -2.0_wp / 3.0_wp ) *&2734 sed_qc(k) = sed_qc_const * nc_sedi**( -2.0_wp / 3.0_wp ) * & 2048 2735 ( qc_1d(k) * hyrho(k) )**( 5.0_wp / 3.0_wp ) * flag 2049 2736 ELSE
Note: See TracChangeset
for help on using the changeset viewer.
