Changeset 1682 for palm/trunk/SOURCE/lpm_collision_kernels.f90
- Timestamp:
- Oct 7, 2015 11:56:08 PM (9 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/lpm_collision_kernels.f90 (modified) (20 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/lpm_collision_kernels.f90
r1520 r1682 1 MODULE lpm_collision_kernels_mod 2 1 !> @file lpm_collision_kernels.f90 3 2 !--------------------------------------------------------------------------------! 4 3 ! This file is part of PALM. … … 20 19 ! Current revisions: 21 20 ! ----------------- 22 ! 21 ! Code annotations made doxygen readable 23 22 ! 24 23 ! Former revisions: … … 89 88 ! Description: 90 89 ! ------------ 91 ! This module calculates collision efficiencies either due to pure gravitational 92 ! effects (Hall kernel, see Hall, 1980: J. Atmos. Sci., 2486-2507) or 93 ! including the effects of (SGS) turbulence (Wang kernel, see Wang and 94 ! Grabowski, 2009: Atmos. Sci. Lett., 10, 1-8). The original code has been 95 ! provided by L.-P. Wang but is substantially reformatted and speed optimized 96 ! here. 97 ! 98 ! ATTENTION: 99 ! Physical quantities (like g, densities, etc.) used in this module still 100 ! have to be adjusted to those values used in the main PALM code. 101 ! Also, quantities in CGS-units should be converted to SI-units eventually. 102 !------------------------------------------------------------------------------! 90 !> This module calculates collision efficiencies either due to pure gravitational 91 !> effects (Hall kernel, see Hall, 1980: J. Atmos. Sci., 2486-2507) or 92 !> including the effects of (SGS) turbulence (Wang kernel, see Wang and 93 !> Grabowski, 2009: Atmos. Sci. Lett., 10, 1-8). The original code has been 94 !> provided by L.-P. Wang but is substantially reformatted and speed optimized 95 !> here. 96 !> 97 !> @attention Physical quantities (like g, densities, etc.) used in this module 98 !> still have to be adjusted to those values used in the main PALM code. 99 !> Also, quantities in CGS-units should be converted to SI-units eventually. 100 !------------------------------------------------------------------------------! 101 MODULE lpm_collision_kernels_mod 102 103 103 104 104 USE constants, & … … 120 120 rclass_lbound, rclass_ubound, recalculate_kernel 121 121 122 REAL(wp) :: epsilon ! :123 REAL(wp) :: eps2 ! :124 REAL(wp) :: rclass_lbound ! :125 REAL(wp) :: rclass_ubound ! :126 REAL(wp) :: urms ! :127 REAL(wp) :: urms2 ! :128 129 REAL(wp), DIMENSION(:), ALLOCATABLE :: epsclass ! :130 REAL(wp), DIMENSION(:), ALLOCATABLE :: radclass ! :131 REAL(wp), DIMENSION(:), ALLOCATABLE :: winf ! :122 REAL(wp) :: epsilon !< 123 REAL(wp) :: eps2 !< 124 REAL(wp) :: rclass_lbound !< 125 REAL(wp) :: rclass_ubound !< 126 REAL(wp) :: urms !< 127 REAL(wp) :: urms2 !< 128 129 REAL(wp), DIMENSION(:), ALLOCATABLE :: epsclass !< 130 REAL(wp), DIMENSION(:), ALLOCATABLE :: radclass !< 131 REAL(wp), DIMENSION(:), ALLOCATABLE :: winf !< 132 132 133 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ec ! :134 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ecf ! :135 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: gck ! :136 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: hkernel ! :137 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: hwratio ! :133 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ec !< 134 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ecf !< 135 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: gck !< 136 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: hkernel !< 137 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: hwratio !< 138 138 139 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ckernel ! :139 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ckernel !< 140 140 141 141 SAVE … … 159 159 160 160 161 !------------------------------------------------------------------------------! 162 ! Description: 163 ! ------------ 164 !> Initialization of the collision efficiency matrix with fixed radius and 165 !> dissipation classes, calculated at simulation start only. 166 !------------------------------------------------------------------------------! 167 161 168 SUBROUTINE init_kernels 162 !------------------------------------------------------------------------------!163 ! Initialization of the collision efficiency matrix with fixed radius and164 ! dissipation classes, calculated at simulation start only.165 !------------------------------------------------------------------------------!166 169 167 170 IMPLICIT NONE 168 171 169 INTEGER(iwp) :: i ! :170 INTEGER(iwp) :: j ! :171 INTEGER(iwp) :: k ! :172 INTEGER(iwp) :: i !< 173 INTEGER(iwp) :: j !< 174 INTEGER(iwp) :: k !< 172 175 173 176 … … 285 288 286 289 !------------------------------------------------------------------------------! 287 ! Calculation of collision kernels during each timestep and for each grid box 290 ! Description: 291 ! ------------ 292 !> Calculation of collision kernels during each timestep and for each grid box 288 293 !------------------------------------------------------------------------------! 289 294 SUBROUTINE recalculate_kernel( i1, j1, k1 ) … … 297 302 IMPLICIT NONE 298 303 299 INTEGER(iwp) :: i ! :300 INTEGER(iwp) :: i1 ! :301 INTEGER(iwp) :: j ! :302 INTEGER(iwp) :: j1 ! :303 INTEGER(iwp) :: k1 ! :304 INTEGER(iwp) :: pend ! :305 INTEGER(iwp) :: pstart ! :304 INTEGER(iwp) :: i !< 305 INTEGER(iwp) :: i1 !< 306 INTEGER(iwp) :: j !< 307 INTEGER(iwp) :: j1 !< 308 INTEGER(iwp) :: k1 !< 309 INTEGER(iwp) :: pend !< 310 INTEGER(iwp) :: pstart !< 306 311 307 312 … … 364 369 365 370 !------------------------------------------------------------------------------! 366 ! Calculation of gck 367 ! This is from Aayala 2008b, page 37ff. 368 ! Necessary input parameters: water density, radii of droplets, air density, 369 ! air viscosity, turbulent dissipation rate, taylor microscale reynolds number, 370 ! gravitational acceleration --> to be replaced by PALM parameters 371 ! Description: 372 ! ------------ 373 !> Calculation of gck 374 !> This is from Aayala 2008b, page 37ff. 375 !> Necessary input parameters: water density, radii of droplets, air density, 376 !> air viscosity, turbulent dissipation rate, taylor microscale reynolds number, 377 !> gravitational acceleration --> to be replaced by PALM parameters 371 378 !------------------------------------------------------------------------------! 372 379 SUBROUTINE turbsd … … 380 387 IMPLICIT NONE 381 388 382 LOGICAL, SAVE :: first = .TRUE. ! :383 384 INTEGER(iwp) :: i ! :385 INTEGER(iwp) :: j ! :386 387 REAL(wp) :: ao ! :388 REAL(wp) :: ao_gr ! :389 REAL(wp) :: bbb ! :390 REAL(wp) :: be ! :391 REAL(wp) :: b1 ! :392 REAL(wp) :: b2 ! :393 REAL(wp) :: ccc ! :394 REAL(wp) :: c1 ! :395 REAL(wp) :: c1_gr ! :396 REAL(wp) :: c2 ! :397 REAL(wp) :: d1 ! :398 REAL(wp) :: d2 ! :399 REAL(wp) :: eta ! :400 REAL(wp) :: e1 ! :401 REAL(wp) :: e2 ! :402 REAL(wp) :: fao_gr ! :403 REAL(wp) :: fr ! :404 REAL(wp) :: grfin ! :405 REAL(wp) :: lambda ! :406 REAL(wp) :: lambda_re ! :407 REAL(wp) :: lf ! :408 REAL(wp) :: rc ! :409 REAL(wp) :: rrp ! :410 REAL(wp) :: sst ! :411 REAL(wp) :: tauk ! :412 REAL(wp) :: tl ! :413 REAL(wp) :: t2 ! :414 REAL(wp) :: tt ! :415 REAL(wp) :: t1 ! :416 REAL(wp) :: vk ! :417 REAL(wp) :: vrms1xy ! :418 REAL(wp) :: vrms2xy ! :419 REAL(wp) :: v1 ! :420 REAL(wp) :: v1v2xy ! :421 REAL(wp) :: v1xysq ! :422 REAL(wp) :: v2 ! :423 REAL(wp) :: v2xysq ! :424 REAL(wp) :: wrfin ! :425 REAL(wp) :: wrgrav2 ! :426 REAL(wp) :: wrtur2xy ! :427 REAL(wp) :: xx ! :428 REAL(wp) :: yy ! :429 REAL(wp) :: z ! :430 431 REAL(wp), DIMENSION(1:radius_classes) :: st ! :432 REAL(wp), DIMENSION(1:radius_classes) :: tau ! :389 LOGICAL, SAVE :: first = .TRUE. !< 390 391 INTEGER(iwp) :: i !< 392 INTEGER(iwp) :: j !< 393 394 REAL(wp) :: ao !< 395 REAL(wp) :: ao_gr !< 396 REAL(wp) :: bbb !< 397 REAL(wp) :: be !< 398 REAL(wp) :: b1 !< 399 REAL(wp) :: b2 !< 400 REAL(wp) :: ccc !< 401 REAL(wp) :: c1 !< 402 REAL(wp) :: c1_gr !< 403 REAL(wp) :: c2 !< 404 REAL(wp) :: d1 !< 405 REAL(wp) :: d2 !< 406 REAL(wp) :: eta !< 407 REAL(wp) :: e1 !< 408 REAL(wp) :: e2 !< 409 REAL(wp) :: fao_gr !< 410 REAL(wp) :: fr !< 411 REAL(wp) :: grfin !< 412 REAL(wp) :: lambda !< 413 REAL(wp) :: lambda_re !< 414 REAL(wp) :: lf !< 415 REAL(wp) :: rc !< 416 REAL(wp) :: rrp !< 417 REAL(wp) :: sst !< 418 REAL(wp) :: tauk !< 419 REAL(wp) :: tl !< 420 REAL(wp) :: t2 !< 421 REAL(wp) :: tt !< 422 REAL(wp) :: t1 !< 423 REAL(wp) :: vk !< 424 REAL(wp) :: vrms1xy !< 425 REAL(wp) :: vrms2xy !< 426 REAL(wp) :: v1 !< 427 REAL(wp) :: v1v2xy !< 428 REAL(wp) :: v1xysq !< 429 REAL(wp) :: v2 !< 430 REAL(wp) :: v2xysq !< 431 REAL(wp) :: wrfin !< 432 REAL(wp) :: wrgrav2 !< 433 REAL(wp) :: wrtur2xy !< 434 REAL(wp) :: xx !< 435 REAL(wp) :: yy !< 436 REAL(wp) :: z !< 437 438 REAL(wp), DIMENSION(1:radius_classes) :: st !< 439 REAL(wp), DIMENSION(1:radius_classes) :: tau !< 433 440 434 441 ! … … 548 555 549 556 !------------------------------------------------------------------------------! 550 ! phi_w as a function 557 ! Description: 558 ! ------------ 559 !> phi_w as a function 551 560 !------------------------------------------------------------------------------! 552 561 REAL(wp) FUNCTION phi_w( a, b, vsett, tau0 ) … … 554 563 IMPLICIT NONE 555 564 556 REAL(wp) :: a ! :557 REAL(wp) :: aa1 ! :558 REAL(wp) :: b ! :559 REAL(wp) :: tau0 ! :560 REAL(wp) :: vsett ! :565 REAL(wp) :: a !< 566 REAL(wp) :: aa1 !< 567 REAL(wp) :: b !< 568 REAL(wp) :: tau0 !< 569 REAL(wp) :: vsett !< 561 570 562 571 aa1 = 1.0_wp / tau0 + 1.0_wp / a + vsett / b … … 567 576 568 577 !------------------------------------------------------------------------------! 569 ! zhi as a function 578 ! Description: 579 ! ------------ 580 !> zhi as a function 570 581 !------------------------------------------------------------------------------! 571 582 REAL(wp) FUNCTION zhi( a, b, vsett1, tau1, vsett2, tau2 ) … … 573 584 IMPLICIT NONE 574 585 575 REAL(wp) :: a ! :576 REAL(wp) :: aa1 ! :577 REAL(wp) :: aa2 ! :578 REAL(wp) :: aa3 ! :579 REAL(wp) :: aa4 ! :580 REAL(wp) :: aa5 ! :581 REAL(wp) :: aa6 ! :582 REAL(wp) :: b ! :583 REAL(wp) :: tau1 ! :584 REAL(wp) :: tau2 ! :585 REAL(wp) :: vsett1 ! :586 REAL(wp) :: vsett2 ! :586 REAL(wp) :: a !< 587 REAL(wp) :: aa1 !< 588 REAL(wp) :: aa2 !< 589 REAL(wp) :: aa3 !< 590 REAL(wp) :: aa4 !< 591 REAL(wp) :: aa5 !< 592 REAL(wp) :: aa6 !< 593 REAL(wp) :: b !< 594 REAL(wp) :: tau1 !< 595 REAL(wp) :: tau2 !< 596 REAL(wp) :: vsett1 !< 597 REAL(wp) :: vsett2 !< 587 598 588 599 aa1 = vsett2 / b - 1.0_wp / tau2 - 1.0_wp / a … … 602 613 603 614 !------------------------------------------------------------------------------! 604 ! Calculation of terminal velocity winf following Equations 10-138 to 10-145 605 ! from (Pruppacher and Klett, 1997) 615 ! Description: 616 ! ------------ 617 !> Calculation of terminal velocity winf following Equations 10-138 to 10-145 618 !> from (Pruppacher and Klett, 1997) 606 619 !------------------------------------------------------------------------------! 607 620 SUBROUTINE fallg … … 619 632 IMPLICIT NONE 620 633 621 INTEGER(iwp) :: i ! :622 INTEGER(iwp) :: j ! :623 624 LOGICAL, SAVE :: first = .TRUE. ! :625 626 REAL(wp), SAVE :: cunh ! :627 REAL(wp), SAVE :: eta ! :628 REAL(wp), SAVE :: phy ! :629 REAL(wp), SAVE :: py ! :630 REAL(wp), SAVE :: rho_a ! :631 REAL(wp), SAVE :: sigma ! :632 REAL(wp), SAVE :: stb ! :633 REAL(wp), SAVE :: stok ! :634 REAL(wp), SAVE :: xlamb ! :635 636 REAL(wp) :: bond ! :637 REAL(wp) :: x ! :638 REAL(wp) :: xrey ! :639 REAL(wp) :: y ! :640 641 REAL(wp), DIMENSION(1:7), SAVE :: b ! :642 REAL(wp), DIMENSION(1:6), SAVE :: c ! :634 INTEGER(iwp) :: i !< 635 INTEGER(iwp) :: j !< 636 637 LOGICAL, SAVE :: first = .TRUE. !< 638 639 REAL(wp), SAVE :: cunh !< 640 REAL(wp), SAVE :: eta !< 641 REAL(wp), SAVE :: phy !< 642 REAL(wp), SAVE :: py !< 643 REAL(wp), SAVE :: rho_a !< 644 REAL(wp), SAVE :: sigma !< 645 REAL(wp), SAVE :: stb !< 646 REAL(wp), SAVE :: stok !< 647 REAL(wp), SAVE :: xlamb !< 648 649 REAL(wp) :: bond !< 650 REAL(wp) :: x !< 651 REAL(wp) :: xrey !< 652 REAL(wp) :: y !< 653 654 REAL(wp), DIMENSION(1:7), SAVE :: b !< 655 REAL(wp), DIMENSION(1:6), SAVE :: c !< 643 656 644 657 ! … … 719 732 720 733 !------------------------------------------------------------------------------! 721 ! Calculation of collision efficiencies for the Hall kernel 734 ! Description: 735 ! ------------ 736 !> Calculation of collision efficiencies for the Hall kernel 722 737 !------------------------------------------------------------------------------! 723 738 SUBROUTINE effic … … 728 743 IMPLICIT NONE 729 744 730 INTEGER(iwp) :: i ! :731 INTEGER(iwp) :: iq ! :732 INTEGER(iwp) :: ir ! :733 INTEGER(iwp) :: j ! :734 INTEGER(iwp) :: k ! :735 736 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ira ! :737 738 LOGICAL, SAVE :: first = .TRUE. ! :739 740 REAL(wp) :: ek ! :741 REAL(wp) :: particle_radius ! :742 REAL(wp) :: pp ! :743 REAL(wp) :: qq ! :744 REAL(wp) :: rq ! :745 746 REAL(wp), DIMENSION(1:21), SAVE :: rat ! :745 INTEGER(iwp) :: i !< 746 INTEGER(iwp) :: iq !< 747 INTEGER(iwp) :: ir !< 748 INTEGER(iwp) :: j !< 749 INTEGER(iwp) :: k !< 750 751 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ira !< 752 753 LOGICAL, SAVE :: first = .TRUE. !< 754 755 REAL(wp) :: ek !< 756 REAL(wp) :: particle_radius !< 757 REAL(wp) :: pp !< 758 REAL(wp) :: qq !< 759 REAL(wp) :: rq !< 760 761 REAL(wp), DIMENSION(1:21), SAVE :: rat !< 747 762 748 REAL(wp), DIMENSION(1:15), SAVE :: r0 ! :763 REAL(wp), DIMENSION(1:15), SAVE :: r0 !< 749 764 750 REAL(wp), DIMENSION(1:15,1:21), SAVE :: ecoll ! :765 REAL(wp), DIMENSION(1:15,1:21), SAVE :: ecoll !< 751 766 752 767 ! … … 888 903 889 904 !------------------------------------------------------------------------------! 890 ! Calculation of enhancement factor for collision efficencies due to turbulence 905 ! Description: 906 ! ------------ 907 !> Calculation of enhancement factor for collision efficencies due to turbulence 891 908 !------------------------------------------------------------------------------! 892 909 SUBROUTINE turb_enhance_eff … … 897 914 IMPLICIT NONE 898 915 899 INTEGER(iwp) :: i ! :900 INTEGER(iwp) :: iq ! :901 INTEGER(iwp) :: ir ! :902 INTEGER(iwp) :: j ! :903 INTEGER(iwp) :: k ! :904 INTEGER(iwp) :: kk ! :905 906 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ira ! :916 INTEGER(iwp) :: i !< 917 INTEGER(iwp) :: iq !< 918 INTEGER(iwp) :: ir !< 919 INTEGER(iwp) :: j !< 920 INTEGER(iwp) :: k !< 921 INTEGER(iwp) :: kk !< 922 923 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: ira !< 907 924 908 LOGICAL, SAVE :: first = .TRUE. ! :909 910 REAL(wp) :: particle_radius ! :911 REAL(wp) :: pp ! :912 REAL(wp) :: qq ! :913 REAL(wp) :: rq ! :914 REAL(wp) :: y1 ! :915 REAL(wp) :: y2 ! :916 REAL(wp) :: y3 ! :917 918 REAL(wp), DIMENSION(1:11), SAVE :: rat ! :925 LOGICAL, SAVE :: first = .TRUE. !< 926 927 REAL(wp) :: particle_radius !< 928 REAL(wp) :: pp !< 929 REAL(wp) :: qq !< 930 REAL(wp) :: rq !< 931 REAL(wp) :: y1 !< 932 REAL(wp) :: y2 !< 933 REAL(wp) :: y3 !< 934 935 REAL(wp), DIMENSION(1:11), SAVE :: rat !< 919 936 920 REAL(wp), DIMENSION(1:7), SAVE :: r0 ! :937 REAL(wp), DIMENSION(1:7), SAVE :: r0 !< 921 938 922 REAL(wp), DIMENSION(1:7,1:11), SAVE :: ecoll_100 ! :923 REAL(wp), DIMENSION(1:7,1:11), SAVE :: ecoll_400 ! :939 REAL(wp), DIMENSION(1:7,1:11), SAVE :: ecoll_100 !< 940 REAL(wp), DIMENSION(1:7,1:11), SAVE :: ecoll_400 !< 924 941 925 942 ! … … 1065 1082 1066 1083 1084 !------------------------------------------------------------------------------! 1085 ! Description: 1086 ! ------------ 1087 !> Collision efficiencies from table 8.2 in Rogers and Yau (1989, 3rd edition). 1088 !> Values are calculated from table by bilinear interpolation. 1089 !------------------------------------------------------------------------------! 1090 1067 1091 SUBROUTINE collision_efficiency_rogers( mean_r, r, e) 1068 !------------------------------------------------------------------------------!1069 ! Collision efficiencies from table 8.2 in Rogers and Yau (1989, 3rd edition).1070 ! Values are calculated from table by bilinear interpolation.1071 !------------------------------------------------------------------------------!1072 1092 1073 1093 IMPLICIT NONE 1074 1094 1075 INTEGER(iwp) :: i ! :1076 INTEGER(iwp) :: j ! :1077 INTEGER(iwp) :: k ! :1078 1079 LOGICAL, SAVE :: first = .TRUE. ! :1080 1081 REAL(wp) :: aa ! :1082 REAL(wp) :: bb ! :1083 REAL(wp) :: cc ! :1084 REAL(wp) :: dd ! :1085 REAL(wp) :: dx ! :1086 REAL(wp) :: dy ! :1087 REAL(wp) :: e ! :1088 REAL(wp) :: gg ! :1089 REAL(wp) :: mean_r ! :1090 REAL(wp) :: mean_rm ! :1091 REAL(wp) :: r ! :1092 REAL(wp) :: rm ! :1093 REAL(wp) :: x ! :1094 REAL(wp) :: y ! :1095 INTEGER(iwp) :: i !< 1096 INTEGER(iwp) :: j !< 1097 INTEGER(iwp) :: k !< 1098 1099 LOGICAL, SAVE :: first = .TRUE. !< 1100 1101 REAL(wp) :: aa !< 1102 REAL(wp) :: bb !< 1103 REAL(wp) :: cc !< 1104 REAL(wp) :: dd !< 1105 REAL(wp) :: dx !< 1106 REAL(wp) :: dy !< 1107 REAL(wp) :: e !< 1108 REAL(wp) :: gg !< 1109 REAL(wp) :: mean_r !< 1110 REAL(wp) :: mean_rm !< 1111 REAL(wp) :: r !< 1112 REAL(wp) :: rm !< 1113 REAL(wp) :: x !< 1114 REAL(wp) :: y !< 1095 1115 1096 REAL(wp), DIMENSION(1:9), SAVE :: collected_r = 0.0_wp ! :1116 REAL(wp), DIMENSION(1:9), SAVE :: collected_r = 0.0_wp !< 1097 1117 1098 REAL(wp), DIMENSION(1:19), SAVE :: collector_r = 0.0_wp ! :1118 REAL(wp), DIMENSION(1:19), SAVE :: collector_r = 0.0_wp !< 1099 1119 1100 REAL(wp), DIMENSION(1:9,1:19), SAVE :: ef = 0.0_wp ! :1120 REAL(wp), DIMENSION(1:9,1:19), SAVE :: ef = 0.0_wp !< 1101 1121 1102 1122 mean_rm = mean_r * 1.0E06_wp
Note: See TracChangeset
for help on using the changeset viewer.
