Changeset 3869 for palm/trunk/SOURCE/bulk_cloud_model_mod.f90
- Timestamp:
- Apr 8, 2019 11:54:20 AM (6 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/bulk_cloud_model_mod.f90 (modified) (20 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/bulk_cloud_model_mod.f90
r3786 r3869 25 25 ! ----------------- 26 26 ! $Id$ 27 ! moving the furniture around ;-) 28 ! 29 ! 3786 2019-03-06 16:58:03Z raasch 27 30 ! unsed variables removed 28 31 ! … … 325 328 bcm_check_data_output, & 326 329 bcm_check_data_output_pr, & 327 bcm_header, &328 330 bcm_init_arrays, & 329 331 bcm_init, & 332 bcm_header, & 333 bcm_actions, & 330 334 bcm_3d_data_averaging, & 331 335 bcm_data_output_2d, & … … 336 340 bcm_wrd_global, & 337 341 bcm_wrd_local, & 338 bcm_actions, &339 342 calc_liquid_water_content 340 343 … … 370 373 END INTERFACE bcm_check_data_output_pr 371 374 375 INTERFACE bcm_init_arrays 376 MODULE PROCEDURE bcm_init_arrays 377 END INTERFACE bcm_init_arrays 378 379 INTERFACE bcm_init 380 MODULE PROCEDURE bcm_init 381 END INTERFACE bcm_init 382 372 383 INTERFACE bcm_header 373 384 MODULE PROCEDURE bcm_header 374 385 END INTERFACE bcm_header 375 376 INTERFACE bcm_init_arrays377 MODULE PROCEDURE bcm_init_arrays378 END INTERFACE bcm_init_arrays379 380 INTERFACE bcm_init381 MODULE PROCEDURE bcm_init382 END INTERFACE bcm_init383 384 INTERFACE bcm_3d_data_averaging385 MODULE PROCEDURE bcm_3d_data_averaging386 END INTERFACE bcm_3d_data_averaging387 388 INTERFACE bcm_data_output_2d389 MODULE PROCEDURE bcm_data_output_2d390 END INTERFACE bcm_data_output_2d391 392 INTERFACE bcm_data_output_3d393 MODULE PROCEDURE bcm_data_output_3d394 END INTERFACE bcm_data_output_3d395 396 INTERFACE bcm_swap_timelevel397 MODULE PROCEDURE bcm_swap_timelevel398 END INTERFACE bcm_swap_timelevel399 400 INTERFACE bcm_rrd_global401 MODULE PROCEDURE bcm_rrd_global402 END INTERFACE bcm_rrd_global403 404 INTERFACE bcm_rrd_local405 MODULE PROCEDURE bcm_rrd_local406 END INTERFACE bcm_rrd_local407 408 INTERFACE bcm_wrd_global409 MODULE PROCEDURE bcm_wrd_global410 END INTERFACE bcm_wrd_global411 412 INTERFACE bcm_wrd_local413 MODULE PROCEDURE bcm_wrd_local414 END INTERFACE bcm_wrd_local415 386 416 387 INTERFACE bcm_actions … … 419 390 END INTERFACE bcm_actions 420 391 421 INTERFACE adjust_cloud 422 MODULE PROCEDURE adjust_cloud 423 MODULE PROCEDURE adjust_cloud_ij 424 END INTERFACE adjust_cloud 425 426 INTERFACE activation 427 MODULE PROCEDURE activation 428 MODULE PROCEDURE activation_ij 429 END INTERFACE activation 430 431 INTERFACE condensation 432 MODULE PROCEDURE condensation 433 MODULE PROCEDURE condensation_ij 434 END INTERFACE condensation 435 436 INTERFACE autoconversion 437 MODULE PROCEDURE autoconversion 438 MODULE PROCEDURE autoconversion_ij 439 END INTERFACE autoconversion 440 441 INTERFACE autoconversion_kessler 442 MODULE PROCEDURE autoconversion_kessler 443 MODULE PROCEDURE autoconversion_kessler_ij 444 END INTERFACE autoconversion_kessler 445 446 INTERFACE accretion 447 MODULE PROCEDURE accretion 448 MODULE PROCEDURE accretion_ij 449 END INTERFACE accretion 450 451 INTERFACE selfcollection_breakup 452 MODULE PROCEDURE selfcollection_breakup 453 MODULE PROCEDURE selfcollection_breakup_ij 454 END INTERFACE selfcollection_breakup 455 456 INTERFACE evaporation_rain 457 MODULE PROCEDURE evaporation_rain 458 MODULE PROCEDURE evaporation_rain_ij 459 END INTERFACE evaporation_rain 460 461 INTERFACE sedimentation_cloud 462 MODULE PROCEDURE sedimentation_cloud 463 MODULE PROCEDURE sedimentation_cloud_ij 464 END INTERFACE sedimentation_cloud 465 466 INTERFACE sedimentation_rain 467 MODULE PROCEDURE sedimentation_rain 468 MODULE PROCEDURE sedimentation_rain_ij 469 END INTERFACE sedimentation_rain 470 471 INTERFACE calc_precipitation_amount 472 MODULE PROCEDURE calc_precipitation_amount 473 MODULE PROCEDURE calc_precipitation_amount_ij 474 END INTERFACE calc_precipitation_amount 475 476 INTERFACE supersaturation 477 MODULE PROCEDURE supersaturation 478 END INTERFACE supersaturation 392 INTERFACE bcm_swap_timelevel 393 MODULE PROCEDURE bcm_swap_timelevel 394 END INTERFACE bcm_swap_timelevel 395 396 INTERFACE bcm_3d_data_averaging 397 MODULE PROCEDURE bcm_3d_data_averaging 398 END INTERFACE bcm_3d_data_averaging 399 400 INTERFACE bcm_data_output_2d 401 MODULE PROCEDURE bcm_data_output_2d 402 END INTERFACE bcm_data_output_2d 403 404 INTERFACE bcm_data_output_3d 405 MODULE PROCEDURE bcm_data_output_3d 406 END INTERFACE bcm_data_output_3d 407 408 INTERFACE bcm_rrd_global 409 MODULE PROCEDURE bcm_rrd_global 410 END INTERFACE bcm_rrd_global 411 412 INTERFACE bcm_rrd_local 413 MODULE PROCEDURE bcm_rrd_local 414 END INTERFACE bcm_rrd_local 415 416 INTERFACE bcm_wrd_global 417 MODULE PROCEDURE bcm_wrd_global 418 END INTERFACE bcm_wrd_global 419 420 INTERFACE bcm_wrd_local 421 MODULE PROCEDURE bcm_wrd_local 422 END INTERFACE bcm_wrd_local 479 423 480 424 INTERFACE calc_liquid_water_content … … 989 933 ! Description: 990 934 ! ------------ 935 !> Header output for bulk cloud module 936 !------------------------------------------------------------------------------! 937 SUBROUTINE bcm_header ( io ) 938 939 940 IMPLICIT NONE 941 942 INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file 943 944 ! 945 !-- Write bulk cloud module header 946 WRITE ( io, 1 ) 947 948 WRITE ( io, 2 ) 949 WRITE ( io, 3 ) 950 951 IF ( microphysics_kessler ) THEN 952 WRITE ( io, 4 ) 'Kessler-Scheme' 953 ENDIF 954 955 IF ( microphysics_seifert ) THEN 956 WRITE ( io, 4 ) 'Seifert-Beheng-Scheme' 957 IF ( cloud_water_sedimentation ) WRITE ( io, 5 ) 958 IF ( collision_turbulence ) WRITE ( io, 6 ) 959 IF ( ventilation_effect ) WRITE ( io, 7 ) 960 IF ( limiter_sedimentation ) WRITE ( io, 8 ) 961 ENDIF 962 963 WRITE ( io, 20 ) 964 WRITE ( io, 21 ) surface_pressure 965 WRITE ( io, 22 ) r_d 966 WRITE ( io, 23 ) rho_surface 967 WRITE ( io, 24 ) c_p 968 WRITE ( io, 25 ) l_v 969 970 IF ( microphysics_seifert ) THEN 971 WRITE ( io, 26 ) 1.0E-6_wp * nc_const 972 WRITE ( io, 27 ) c_sedimentation 973 ENDIF 974 975 976 1 FORMAT ( //' Bulk cloud module information:'/ & 977 ' ------------------------------------------'/ ) 978 2 FORMAT ( '--> Bulk scheme with liquid water potential temperature and'/ & 979 ' total water content is used.' ) 980 3 FORMAT ( '--> Condensation is parameterized via 0% - or 100% scheme.' ) 981 4 FORMAT ( '--> Precipitation parameterization via ', A ) 982 983 5 FORMAT ( '--> Cloud water sedimentation parameterization via Stokes law' ) 984 6 FORMAT ( '--> Turbulence effects on precipitation process' ) 985 7 FORMAT ( '--> Ventilation effects on evaporation of rain drops' ) 986 8 FORMAT ( '--> Slope limiter used for sedimentation process' ) 987 988 20 FORMAT ( '--> Essential parameters:' ) 989 21 FORMAT ( ' Surface pressure : p_0 = ', F7.2, ' hPa') 990 22 FORMAT ( ' Gas constant : R = ', F5.1, ' J/(kg K)') 991 23 FORMAT ( ' Density of air : rho_0 = ', F6.3, ' kg/m**3') 992 24 FORMAT ( ' Specific heat cap. : c_p = ', F6.1, ' J/(kg K)') 993 25 FORMAT ( ' Vapourization heat : L_v = ', E9.2, ' J/kg') 994 26 FORMAT ( ' Droplet density : N_c = ', F6.1, ' 1/cm**3' ) 995 27 FORMAT ( ' Sedimentation Courant number : C_s = ', F4.1 ) 996 997 998 END SUBROUTINE bcm_header 999 1000 1001 !------------------------------------------------------------------------------! 1002 ! Description: 1003 ! ------------ 1004 !> Control of microphysics for all grid points 1005 !------------------------------------------------------------------------------! 1006 SUBROUTINE bcm_actions 1007 1008 IMPLICIT NONE 1009 1010 IF ( large_scale_forcing .AND. lsf_surf ) THEN 1011 ! 1012 !-- Calculate vertical profile of the hydrostatic pressure (hyp) 1013 hyp = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp) 1014 d_exner = exner_function_invers(hyp) 1015 exner = 1.0_wp / exner_function_invers(hyp) 1016 hyrho = ideal_gas_law_rho_pt(hyp, pt_init) 1017 ! 1018 !-- Compute reference density 1019 rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp)) 1020 ENDIF 1021 1022 ! 1023 !-- Compute length of time step 1024 IF ( call_microphysics_at_all_substeps ) THEN 1025 dt_micro = dt_3d * weight_pres(intermediate_timestep_count) 1026 ELSE 1027 dt_micro = dt_3d 1028 ENDIF 1029 1030 ! 1031 !-- Reset precipitation rate 1032 IF ( intermediate_timestep_count == 1 ) prr = 0.0_wp 1033 1034 ! 1035 !-- Compute cloud physics 1036 IF ( microphysics_kessler ) THEN 1037 1038 CALL autoconversion_kessler 1039 IF ( cloud_water_sedimentation ) CALL sedimentation_cloud 1040 1041 ELSEIF ( microphysics_seifert ) THEN 1042 1043 CALL adjust_cloud 1044 IF ( microphysics_morrison ) CALL activation 1045 IF ( microphysics_morrison ) CALL condensation 1046 CALL autoconversion 1047 CALL accretion 1048 CALL selfcollection_breakup 1049 CALL evaporation_rain 1050 CALL sedimentation_rain 1051 IF ( cloud_water_sedimentation ) CALL sedimentation_cloud 1052 1053 ENDIF 1054 1055 CALL calc_precipitation_amount 1056 1057 END SUBROUTINE bcm_actions 1058 1059 1060 !------------------------------------------------------------------------------! 1061 ! Description: 1062 ! ------------ 1063 !> Control of microphysics for grid points i,j 1064 !------------------------------------------------------------------------------! 1065 1066 SUBROUTINE bcm_actions_ij( i, j ) 1067 1068 IMPLICIT NONE 1069 1070 INTEGER(iwp) :: i !< 1071 INTEGER(iwp) :: j !< 1072 1073 IF ( large_scale_forcing .AND. lsf_surf ) THEN 1074 ! 1075 !-- Calculate vertical profile of the hydrostatic pressure (hyp) 1076 hyp = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp) 1077 d_exner = exner_function_invers(hyp) 1078 exner = 1.0_wp / exner_function_invers(hyp) 1079 hyrho = ideal_gas_law_rho_pt(hyp, pt_init) 1080 ! 1081 !-- Compute reference density 1082 rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp)) 1083 ENDIF 1084 1085 ! 1086 !-- Compute length of time step 1087 IF ( call_microphysics_at_all_substeps ) THEN 1088 dt_micro = dt_3d * weight_pres(intermediate_timestep_count) 1089 ELSE 1090 dt_micro = dt_3d 1091 ENDIF 1092 ! 1093 !-- Reset precipitation rate 1094 IF ( intermediate_timestep_count == 1 ) prr(:,j,i) = 0.0_wp 1095 1096 ! 1097 !-- Compute cloud physics 1098 IF( microphysics_kessler ) THEN 1099 1100 CALL autoconversion_kessler_ij( i,j ) 1101 IF ( cloud_water_sedimentation ) CALL sedimentation_cloud_ij( i,j ) 1102 1103 ELSEIF ( microphysics_seifert ) THEN 1104 1105 CALL adjust_cloud_ij( i,j ) 1106 IF ( microphysics_morrison ) CALL activation_ij( i,j ) 1107 IF ( microphysics_morrison ) CALL condensation_ij( i,j ) 1108 CALL autoconversion_ij( i,j ) 1109 CALL accretion_ij( i,j ) 1110 CALL selfcollection_breakup_ij( i,j ) 1111 CALL evaporation_rain_ij( i,j ) 1112 CALL sedimentation_rain_ij( i,j ) 1113 IF ( cloud_water_sedimentation ) CALL sedimentation_cloud_ij( i,j ) 1114 1115 ENDIF 1116 1117 CALL calc_precipitation_amount_ij( i,j ) 1118 1119 END SUBROUTINE bcm_actions_ij 1120 1121 1122 !------------------------------------------------------------------------------! 1123 ! Description: 1124 ! ------------ 991 1125 !> Swapping of timelevels 992 1126 !------------------------------------------------------------------------------! … … 1028 1162 1029 1163 END SUBROUTINE bcm_swap_timelevel 1030 1031 1032 !------------------------------------------------------------------------------!1033 ! Description:1034 ! ------------1035 !> Header output for bulk cloud module1036 !------------------------------------------------------------------------------!1037 SUBROUTINE bcm_header ( io )1038 1039 1040 IMPLICIT NONE1041 1042 INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file1043 1044 !1045 !-- Write bulk cloud module header1046 WRITE ( io, 1 )1047 1048 WRITE ( io, 2 )1049 WRITE ( io, 3 )1050 1051 IF ( microphysics_kessler ) THEN1052 WRITE ( io, 4 ) 'Kessler-Scheme'1053 ENDIF1054 1055 IF ( microphysics_seifert ) THEN1056 WRITE ( io, 4 ) 'Seifert-Beheng-Scheme'1057 IF ( cloud_water_sedimentation ) WRITE ( io, 5 )1058 IF ( collision_turbulence ) WRITE ( io, 6 )1059 IF ( ventilation_effect ) WRITE ( io, 7 )1060 IF ( limiter_sedimentation ) WRITE ( io, 8 )1061 ENDIF1062 1063 WRITE ( io, 20 )1064 WRITE ( io, 21 ) surface_pressure1065 WRITE ( io, 22 ) r_d1066 WRITE ( io, 23 ) rho_surface1067 WRITE ( io, 24 ) c_p1068 WRITE ( io, 25 ) l_v1069 1070 IF ( microphysics_seifert ) THEN1071 WRITE ( io, 26 ) 1.0E-6_wp * nc_const1072 WRITE ( io, 27 ) c_sedimentation1073 ENDIF1074 1075 1076 1 FORMAT ( //' Bulk cloud module information:'/ &1077 ' ------------------------------------------'/ )1078 2 FORMAT ( '--> Bulk scheme with liquid water potential temperature and'/ &1079 ' total water content is used.' )1080 3 FORMAT ( '--> Condensation is parameterized via 0% - or 100% scheme.' )1081 4 FORMAT ( '--> Precipitation parameterization via ', A )1082 1083 5 FORMAT ( '--> Cloud water sedimentation parameterization via Stokes law' )1084 6 FORMAT ( '--> Turbulence effects on precipitation process' )1085 7 FORMAT ( '--> Ventilation effects on evaporation of rain drops' )1086 8 FORMAT ( '--> Slope limiter used for sedimentation process' )1087 1088 20 FORMAT ( '--> Essential parameters:' )1089 21 FORMAT ( ' Surface pressure : p_0 = ', F7.2, ' hPa')1090 22 FORMAT ( ' Gas constant : R = ', F5.1, ' J/(kg K)')1091 23 FORMAT ( ' Density of air : rho_0 = ', F6.3, ' kg/m**3')1092 24 FORMAT ( ' Specific heat cap. : c_p = ', F6.1, ' J/(kg K)')1093 25 FORMAT ( ' Vapourization heat : L_v = ', E9.2, ' J/kg')1094 26 FORMAT ( ' Droplet density : N_c = ', F6.1, ' 1/cm**3' )1095 27 FORMAT ( ' Sedimentation Courant number : C_s = ', F4.1 )1096 1097 1098 END SUBROUTINE bcm_header1099 1164 1100 1165 … … 1993 2058 END SUBROUTINE bcm_wrd_local 1994 2059 1995 1996 !------------------------------------------------------------------------------!1997 ! Description:1998 ! ------------1999 !> Control of microphysics for all grid points2000 !------------------------------------------------------------------------------!2001 SUBROUTINE bcm_actions2002 2003 IMPLICIT NONE2004 2005 IF ( large_scale_forcing .AND. lsf_surf ) THEN2006 !2007 !-- Calculate vertical profile of the hydrostatic pressure (hyp)2008 hyp = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)2009 d_exner = exner_function_invers(hyp)2010 exner = 1.0_wp / exner_function_invers(hyp)2011 hyrho = ideal_gas_law_rho_pt(hyp, pt_init)2012 !2013 !-- Compute reference density2014 rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp))2015 ENDIF2016 2017 !2018 !-- Compute length of time step2019 IF ( call_microphysics_at_all_substeps ) THEN2020 dt_micro = dt_3d * weight_pres(intermediate_timestep_count)2021 ELSE2022 dt_micro = dt_3d2023 ENDIF2024 2025 !2026 !-- Reset precipitation rate2027 IF ( intermediate_timestep_count == 1 ) prr = 0.0_wp2028 2029 !2030 !-- Compute cloud physics2031 IF ( microphysics_kessler ) THEN2032 2033 CALL autoconversion_kessler2034 IF ( cloud_water_sedimentation ) CALL sedimentation_cloud2035 2036 ELSEIF ( microphysics_seifert ) THEN2037 2038 CALL adjust_cloud2039 IF ( microphysics_morrison ) CALL activation2040 IF ( microphysics_morrison ) CALL condensation2041 CALL autoconversion2042 CALL accretion2043 CALL selfcollection_breakup2044 CALL evaporation_rain2045 CALL sedimentation_rain2046 IF ( cloud_water_sedimentation ) CALL sedimentation_cloud2047 2048 ENDIF2049 2050 CALL calc_precipitation_amount2051 2052 END SUBROUTINE bcm_actions2053 2054 2060 !------------------------------------------------------------------------------! 2055 2061 ! Description: … … 2107 2113 2108 2114 END SUBROUTINE adjust_cloud 2115 2116 !------------------------------------------------------------------------------! 2117 ! Description: 2118 ! ------------ 2119 !> Adjust number of raindrops to avoid nonlinear effects in 2120 !> sedimentation and evaporation of rain drops due to too small or 2121 !> too big weights of rain drops (Stevens and Seifert, 2008). 2122 !> The same procedure is applied to cloud droplets if they are determined 2123 !> prognostically. Call for grid point i,j 2124 !------------------------------------------------------------------------------! 2125 SUBROUTINE adjust_cloud_ij( i, j ) 2126 2127 IMPLICIT NONE 2128 2129 INTEGER(iwp) :: i !< 2130 INTEGER(iwp) :: j !< 2131 INTEGER(iwp) :: k !< 2132 2133 REAL(wp) :: flag !< flag to indicate first grid level above surface 2134 2135 DO k = nzb+1, nzt 2136 ! 2137 !-- Predetermine flag to mask topography 2138 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 2139 2140 IF ( qr(k,j,i) <= eps_sb ) THEN 2141 qr(k,j,i) = 0.0_wp 2142 nr(k,j,i) = 0.0_wp 2143 ELSE 2144 ! 2145 !-- Adjust number of raindrops to avoid nonlinear effects in 2146 !-- sedimentation and evaporation of rain drops due to too small or 2147 !-- too big weights of rain drops (Stevens and Seifert, 2008). 2148 IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) ) THEN 2149 nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag 2150 ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) ) THEN 2151 nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag 2152 ENDIF 2153 2154 ENDIF 2155 2156 IF ( microphysics_morrison ) THEN 2157 IF ( qc(k,j,i) <= eps_sb ) THEN 2158 qc(k,j,i) = 0.0_wp 2159 nc(k,j,i) = 0.0_wp 2160 ELSE 2161 IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) ) THEN 2162 nc(k,j,i) = qc(k,j,i) * hyrho(k) / xamin * flag 2163 ENDIF 2164 ENDIF 2165 ENDIF 2166 2167 ENDDO 2168 2169 END SUBROUTINE adjust_cloud_ij 2109 2170 2110 2171 !------------------------------------------------------------------------------! … … 2214 2275 END SUBROUTINE activation 2215 2276 2277 !------------------------------------------------------------------------------! 2278 ! Description: 2279 ! ------------ 2280 !> Calculate number of activated condensation nucleii after simple activation 2281 !> scheme of Twomey, 1959. 2282 !------------------------------------------------------------------------------! 2283 SUBROUTINE activation_ij( i, j ) 2284 2285 IMPLICIT NONE 2286 2287 INTEGER(iwp) :: i !< 2288 INTEGER(iwp) :: j !< 2289 INTEGER(iwp) :: k !< 2290 2291 REAL(wp) :: activ !< 2292 REAL(wp) :: afactor !< 2293 REAL(wp) :: beta_act !< 2294 REAL(wp) :: bfactor !< 2295 REAL(wp) :: flag !< flag to indicate first grid level above surface 2296 REAL(wp) :: k_act !< 2297 REAL(wp) :: n_act !< 2298 REAL(wp) :: n_ccn !< 2299 REAL(wp) :: s_0 !< 2300 REAL(wp) :: sat_max !< 2301 REAL(wp) :: sigma !< 2302 REAL(wp) :: sigma_act !< 2303 2304 DO k = nzb+1, nzt 2305 ! 2306 !-- Predetermine flag to mask topography 2307 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 2308 ! 2309 !-- Call calculation of supersaturation 2310 CALL supersaturation ( i, j, k ) 2311 ! 2312 !-- Prescribe parameters for activation 2313 !-- (see: Bott + Trautmann, 2002, Atm. Res., 64) 2314 k_act = 0.7_wp 2315 activ = 0.0_wp 2316 2317 IF ( sat > 0.0 .AND. .NOT. curvature_solution_effects_bulk ) THEN 2318 ! 2319 !-- Compute the number of activated Aerosols 2320 !-- (see: Twomey, 1959, Pure and applied Geophysics, 43) 2321 n_act = na_init * sat**k_act 2322 ! 2323 !-- Compute the number of cloud droplets 2324 !-- (see: Morrison + Grabowski, 2007, JAS, 64) 2325 ! activ = MAX( n_act - nc_d1(k), 0.0_wp) / dt_micro 2326 2327 ! 2328 !-- Compute activation rate after Khairoutdinov and Kogan 2329 !-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128) 2330 sat_max = 0.8_wp / 100.0_wp 2331 activ = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN & 2332 ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) / & 2333 dt_micro 2334 2335 nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro), na_init) 2336 ELSEIF ( sat > 0.0 .AND. curvature_solution_effects_bulk ) THEN 2337 ! 2338 !-- Curvature effect (afactor) with surface tension 2339 !-- parameterization by Straka (2009) 2340 sigma = 0.0761_wp - 0.000155_wp * ( t_l - 273.15_wp ) 2341 afactor = 2.0_wp * sigma / ( rho_l * r_v * t_l ) 2342 ! 2343 !-- Solute effect (bfactor) 2344 bfactor = vanthoff * molecular_weight_of_water * & 2345 rho_s / ( molecular_weight_of_solute * rho_l ) 2346 2347 ! 2348 !-- Prescribe power index that describes the soluble fraction 2349 !-- of an aerosol particle (beta). 2350 !-- (see: Morrison + Grabowski, 2007, JAS, 64) 2351 beta_act = 0.5_wp 2352 sigma_act = sigma_bulk**( 1.0_wp + beta_act ) 2353 ! 2354 !-- Calculate mean geometric supersaturation (s_0) with 2355 !-- parameterization by Khvorostyanov and Curry (2006) 2356 s_0 = dry_aerosol_radius **(- ( 1.0_wp + beta_act ) ) * & 2357 ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp 2358 2359 ! 2360 !-- Calculate number of activated CCN as a function of 2361 !-- supersaturation and taking Koehler theory into account 2362 !-- (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111) 2363 n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF( & 2364 LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) ) 2365 activ = MAX( ( n_ccn ) / dt_micro, 0.0_wp ) 2366 2367 nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro * flag), na_init) 2368 ENDIF 2369 2370 ENDDO 2371 2372 END SUBROUTINE activation_ij 2373 2216 2374 2217 2375 !------------------------------------------------------------------------------! … … 2297 2455 2298 2456 END SUBROUTINE condensation 2457 2458 !------------------------------------------------------------------------------! 2459 ! Description: 2460 ! ------------ 2461 !> Calculate condensation rate for cloud water content (after Khairoutdinov and 2462 !> Kogan, 2000). 2463 !------------------------------------------------------------------------------! 2464 SUBROUTINE condensation_ij( i, j ) 2465 2466 IMPLICIT NONE 2467 2468 INTEGER(iwp) :: i !< 2469 INTEGER(iwp) :: j !< 2470 INTEGER(iwp) :: k !< 2471 2472 REAL(wp) :: cond !< 2473 REAL(wp) :: cond_max !< 2474 REAL(wp) :: dc !< 2475 REAL(wp) :: evap !< 2476 REAL(wp) :: flag !< flag to indicate first grid level above surface 2477 REAL(wp) :: g_fac !< 2478 REAL(wp) :: nc_0 !< 2479 REAL(wp) :: temp !< 2480 REAL(wp) :: xc !< 2481 2482 2483 DO k = nzb+1, nzt 2484 ! 2485 !-- Predetermine flag to mask topography 2486 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 2487 ! 2488 !-- Call calculation of supersaturation 2489 CALL supersaturation ( i, j, k ) 2490 ! 2491 !-- Actual temperature: 2492 temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) ) 2493 2494 g_fac = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * & 2495 l_v / ( thermal_conductivity_l * temp ) & 2496 + r_v * temp / ( diff_coeff_l * e_s ) & 2497 ) 2498 ! 2499 !-- Mean weight of cloud drops 2500 IF ( nc(k,j,i) <= 0.0_wp) CYCLE 2501 xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin) 2502 ! 2503 !-- Weight averaged diameter of cloud drops: 2504 dc = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp ) 2505 ! 2506 !-- Integral diameter of cloud drops 2507 nc_0 = nc(k,j,i) * dc 2508 ! 2509 !-- Condensation needs only to be calculated in supersaturated regions 2510 IF ( sat > 0.0_wp ) THEN 2511 ! 2512 !-- Condensation rate of cloud water content 2513 !-- after KK scheme. 2514 !-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128) 2515 cond = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k) 2516 cond_max = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i) 2517 cond = MIN( cond, cond_max / dt_micro ) 2518 2519 qc(k,j,i) = qc(k,j,i) + cond * dt_micro * flag 2520 ELSEIF ( sat < 0.0_wp ) THEN 2521 evap = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k) 2522 evap = MAX( evap, -qc(k,j,i) / dt_micro ) 2523 2524 qc(k,j,i) = qc(k,j,i) + evap * dt_micro * flag 2525 ENDIF 2526 ENDDO 2527 2528 END SUBROUTINE condensation_ij 2299 2529 2300 2530 … … 2432 2662 ! Description: 2433 2663 ! ------------ 2664 !> Autoconversion rate (Seifert and Beheng, 2006). Call for grid point i,j 2665 !------------------------------------------------------------------------------! 2666 SUBROUTINE autoconversion_ij( i, j ) 2667 2668 IMPLICIT NONE 2669 2670 INTEGER(iwp) :: i !< 2671 INTEGER(iwp) :: j !< 2672 INTEGER(iwp) :: k !< 2673 2674 REAL(wp) :: alpha_cc !< 2675 REAL(wp) :: autocon !< 2676 REAL(wp) :: dissipation !< 2677 REAL(wp) :: flag !< flag to indicate first grid level above surface 2678 REAL(wp) :: k_au !< 2679 REAL(wp) :: l_mix !< 2680 REAL(wp) :: nc_auto !< 2681 REAL(wp) :: nu_c !< 2682 REAL(wp) :: phi_au !< 2683 REAL(wp) :: r_cc !< 2684 REAL(wp) :: rc !< 2685 REAL(wp) :: re_lambda !< 2686 REAL(wp) :: sigma_cc !< 2687 REAL(wp) :: tau_cloud !< 2688 REAL(wp) :: xc !< 2689 2690 DO k = nzb+1, nzt 2691 ! 2692 !-- Predetermine flag to mask topography 2693 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 2694 IF ( microphysics_morrison ) THEN 2695 nc_auto = nc(k,j,i) 2696 ELSE 2697 nc_auto = nc_const 2698 ENDIF 2699 2700 IF ( qc(k,j,i) > eps_sb .AND. nc_auto > eps_mr ) THEN 2701 2702 k_au = k_cc / ( 20.0_wp * x0 ) 2703 ! 2704 !-- Intern time scale of coagulation (Seifert and Beheng, 2006): 2705 !-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )) 2706 tau_cloud = MAX( 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) ), & 2707 0.0_wp ) 2708 ! 2709 !-- Universal function for autoconversion process 2710 !-- (Seifert and Beheng, 2006): 2711 phi_au = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**3 2712 ! 2713 !-- Shape parameter of gamma distribution (Geoffroy et al., 2010): 2714 !-- (Use constant nu_c = 1.0_wp instead?) 2715 nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp ) 2716 ! 2717 !-- Mean weight of cloud droplets: 2718 xc = hyrho(k) * qc(k,j,i) / nc_auto 2719 ! 2720 !-- Parameterized turbulence effects on autoconversion (Seifert, 2721 !-- Nuijens and Stevens, 2010) 2722 IF ( collision_turbulence ) THEN 2723 ! 2724 !-- Weight averaged radius of cloud droplets: 2725 rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp ) 2726 2727 alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c ) 2728 r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c ) 2729 sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c ) 2730 ! 2731 !-- Mixing length (neglecting distance to ground and stratification) 2732 l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp ) 2733 ! 2734 !-- Limit dissipation rate according to Seifert, Nuijens and 2735 !-- Stevens (2010) 2736 dissipation = MIN( 0.06_wp, diss(k,j,i) ) 2737 ! 2738 !-- Compute Taylor-microscale Reynolds number: 2739 re_lambda = 6.0_wp / 11.0_wp * & 2740 ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) * & 2741 SQRT( 15.0_wp / kin_vis_air ) * & 2742 dissipation**( 1.0_wp / 6.0_wp ) 2743 ! 2744 !-- The factor of 1.0E4 is needed to convert the dissipation rate 2745 !-- from m2 s-3 to cm2 s-3. 2746 k_au = k_au * ( 1.0_wp + & 2747 dissipation * 1.0E4_wp * & 2748 ( re_lambda * 1.0E-3_wp )**0.25_wp * & 2749 ( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) / & 2750 sigma_cc )**2 & 2751 ) + beta_cc & 2752 ) & 2753 ) 2754 ENDIF 2755 ! 2756 !-- Autoconversion rate (Seifert and Beheng, 2006): 2757 autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / & 2758 ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 * & 2759 ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) * & 2760 rho_surface 2761 autocon = MIN( autocon, qc(k,j,i) / dt_micro ) 2762 2763 qr(k,j,i) = qr(k,j,i) + autocon * dt_micro * flag 2764 qc(k,j,i) = qc(k,j,i) - autocon * dt_micro * flag 2765 nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro * flag 2766 IF ( microphysics_morrison ) THEN 2767 nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp * & 2768 autocon / x0 * hyrho(k) * dt_micro * flag ) 2769 ENDIF 2770 2771 ENDIF 2772 2773 ENDDO 2774 2775 END SUBROUTINE autoconversion_ij 2776 2777 2778 !------------------------------------------------------------------------------! 2779 ! Description: 2780 ! ------------ 2434 2781 !> Autoconversion process (Kessler, 1969). 2435 2782 !------------------------------------------------------------------------------! … … 2477 2824 2478 2825 END SUBROUTINE autoconversion_kessler 2826 2827 !------------------------------------------------------------------------------! 2828 ! Description: 2829 ! ------------ 2830 !> Autoconversion process (Kessler, 1969). 2831 !------------------------------------------------------------------------------! 2832 SUBROUTINE autoconversion_kessler_ij( i, j ) 2833 2834 2835 IMPLICIT NONE 2836 2837 INTEGER(iwp) :: i !< 2838 INTEGER(iwp) :: j !< 2839 INTEGER(iwp) :: k !< 2840 INTEGER(iwp) :: k_wall !< topography top index 2841 2842 REAL(wp) :: dqdt_precip !< 2843 REAL(wp) :: flag !< flag to indicate first grid level above surface 2844 2845 ! 2846 !-- Determine vertical index of topography top 2847 k_wall = get_topography_top_index_ji( j, i, 's' ) 2848 DO k = nzb+1, nzt 2849 ! 2850 !-- Predetermine flag to mask topography 2851 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 2852 2853 IF ( qc(k,j,i) > ql_crit ) THEN 2854 dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit ) 2855 ELSE 2856 dqdt_precip = 0.0_wp 2857 ENDIF 2858 2859 qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro * flag 2860 q(k,j,i) = q(k,j,i) - dqdt_precip * dt_micro * flag 2861 pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * lv_d_cp * d_exner(k) & 2862 * flag 2863 2864 ! 2865 !-- Compute the rain rate (stored on surface grid point) 2866 prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag 2867 2868 ENDDO 2869 2870 END SUBROUTINE autoconversion_kessler_ij 2479 2871 2480 2872 … … 2564 2956 END SUBROUTINE accretion 2565 2957 2958 !------------------------------------------------------------------------------! 2959 ! Description: 2960 ! ------------ 2961 !> Accretion rate (Seifert and Beheng, 2006). Call for grid point i,j 2962 !------------------------------------------------------------------------------! 2963 SUBROUTINE accretion_ij( i, j ) 2964 2965 IMPLICIT NONE 2966 2967 INTEGER(iwp) :: i !< 2968 INTEGER(iwp) :: j !< 2969 INTEGER(iwp) :: k !< 2970 2971 REAL(wp) :: accr !< 2972 REAL(wp) :: flag !< flag to indicate first grid level above surface 2973 REAL(wp) :: k_cr !< 2974 REAL(wp) :: nc_accr !< 2975 REAL(wp) :: phi_ac !< 2976 REAL(wp) :: tau_cloud !< 2977 REAL(wp) :: xc !< 2978 2979 2980 DO k = nzb+1, nzt 2981 ! 2982 !-- Predetermine flag to mask topography 2983 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 2984 IF ( microphysics_morrison ) THEN 2985 nc_accr = nc(k,j,i) 2986 ELSE 2987 nc_accr = nc_const 2988 ENDIF 2989 2990 IF ( ( qc(k,j,i) > eps_sb ) .AND. ( qr(k,j,i) > eps_sb ) .AND. & 2991 ( nc_accr > eps_mr ) ) THEN 2992 ! 2993 !-- Intern time scale of coagulation (Seifert and Beheng, 2006): 2994 tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ) 2995 ! 2996 !-- Universal function for accretion process 2997 !-- (Seifert and Beheng, 2001): 2998 phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4 2999 3000 ! 3001 !-- Mean weight of cloud drops 3002 xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin) 3003 ! 3004 !-- Parameterized turbulence effects on autoconversion (Seifert, 3005 !-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to 3006 !-- convert the dissipation rate (diss) from m2 s-3 to cm2 s-3. 3007 IF ( collision_turbulence ) THEN 3008 k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * & 3009 MIN( 600.0_wp, & 3010 diss(k,j,i) * 1.0E4_wp )**0.25_wp & 3011 ) 3012 ELSE 3013 k_cr = k_cr0 3014 ENDIF 3015 ! 3016 !-- Accretion rate (Seifert and Beheng, 2006): 3017 accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac * & 3018 SQRT( rho_surface * hyrho(k) ) 3019 accr = MIN( accr, qc(k,j,i) / dt_micro ) 3020 3021 qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag 3022 qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag 3023 IF ( microphysics_morrison ) THEN 3024 nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), accr / xc * & 3025 hyrho(k) * dt_micro * flag & 3026 ) 3027 ENDIF 3028 3029 3030 ENDIF 3031 3032 ENDDO 3033 3034 END SUBROUTINE accretion_ij 3035 2566 3036 2567 3037 !------------------------------------------------------------------------------! … … 2622 3092 2623 3093 END SUBROUTINE selfcollection_breakup 3094 3095 3096 !------------------------------------------------------------------------------! 3097 ! Description: 3098 ! ------------ 3099 !> Collisional breakup rate (Seifert, 2008). Call for grid point i,j 3100 !------------------------------------------------------------------------------! 3101 SUBROUTINE selfcollection_breakup_ij( i, j ) 3102 3103 IMPLICIT NONE 3104 3105 INTEGER(iwp) :: i !< 3106 INTEGER(iwp) :: j !< 3107 INTEGER(iwp) :: k !< 3108 3109 REAL(wp) :: breakup !< 3110 REAL(wp) :: dr !< 3111 REAL(wp) :: flag !< flag to indicate first grid level above surface 3112 REAL(wp) :: phi_br !< 3113 REAL(wp) :: selfcoll !< 3114 3115 DO k = nzb+1, nzt 3116 ! 3117 !-- Predetermine flag to mask topography 3118 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 3119 3120 IF ( qr(k,j,i) > eps_sb ) THEN 3121 ! 3122 !-- Selfcollection rate (Seifert and Beheng, 2001): 3123 selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * SQRT( hyrho(k) * rho_surface ) 3124 ! 3125 !-- Weight averaged diameter of rain drops: 3126 dr = ( hyrho(k) * qr(k,j,i) / nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp ) 3127 ! 3128 !-- Collisional breakup rate (Seifert, 2008): 3129 IF ( dr >= 0.3E-3_wp ) THEN 3130 phi_br = k_br * ( dr - 1.1E-3_wp ) 3131 breakup = selfcoll * ( phi_br + 1.0_wp ) 3132 ELSE 3133 breakup = 0.0_wp 3134 ENDIF 3135 3136 selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro ) 3137 nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro * flag 3138 3139 ENDIF 3140 ENDDO 3141 3142 END SUBROUTINE selfcollection_breakup_ij 2624 3143 2625 3144 … … 2753 3272 ! Description: 2754 3273 ! ------------ 2755 !> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR).2756 !------------------------------------------------------------------------------!2757 SUBROUTINE sedimentation_cloud2758 2759 2760 IMPLICIT NONE2761 2762 INTEGER(iwp) :: i !<2763 INTEGER(iwp) :: j !<2764 INTEGER(iwp) :: k !<2765 2766 REAL(wp) :: flag !< flag to mask topography grid points2767 REAL(wp) :: nc_sedi !<2768 2769 REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !<2770 REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nc !<2771 2772 2773 CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' )2774 2775 sed_qc(nzt+1) = 0.0_wp2776 sed_nc(nzt+1) = 0.0_wp2777 2778 DO i = nxlg, nxrg2779 DO j = nysg, nyng2780 DO k = nzt, nzb+1, -12781 !2782 !-- Predetermine flag to mask topography2783 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )2784 2785 IF ( microphysics_morrison ) THEN2786 nc_sedi = nc(k,j,i)2787 ELSE2788 nc_sedi = nc_const2789 ENDIF2790 2791 !2792 !-- Sedimentation fluxes for number concentration are only calculated2793 !-- for cloud_scheme = 'morrison'2794 IF ( microphysics_morrison ) THEN2795 IF ( qc(k,j,i) > eps_sb .AND. nc(k,j,i) > eps_mr ) THEN2796 sed_nc(k) = sed_qc_const * &2797 ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) * &2798 ( nc(k,j,i) )**( 1.0_wp / 3.0_wp )2799 ELSE2800 sed_nc(k) = 0.0_wp2801 ENDIF2802 2803 sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) * &2804 nc(k,j,i) / dt_micro + sed_nc(k+1) &2805 ) * flag2806 2807 nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) * &2808 ddzu(k+1) / hyrho(k) * dt_micro * flag2809 ENDIF2810 2811 IF ( qc(k,j,i) > eps_sb .AND. nc_sedi > eps_mr ) THEN2812 sed_qc(k) = sed_qc_const * nc_sedi**( -2.0_wp / 3.0_wp ) * &2813 ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp ) * &2814 flag2815 ELSE2816 sed_qc(k) = 0.0_wp2817 ENDIF2818 2819 sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) / &2820 dt_micro + sed_qc(k+1) &2821 ) * flag2822 2823 q(k,j,i) = q(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * &2824 ddzu(k+1) / hyrho(k) * dt_micro * flag2825 qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * &2826 ddzu(k+1) / hyrho(k) * dt_micro * flag2827 pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) * &2828 ddzu(k+1) / hyrho(k) * lv_d_cp * &2829 d_exner(k) * dt_micro * flag2830 2831 !2832 !-- Compute the precipitation rate due to cloud (fog) droplets2833 IF ( call_microphysics_at_all_substeps ) THEN2834 prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) &2835 * weight_substep(intermediate_timestep_count) &2836 * flag2837 ELSE2838 prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag2839 ENDIF2840 2841 ENDDO2842 ENDDO2843 ENDDO2844 2845 CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' )2846 2847 END SUBROUTINE sedimentation_cloud2848 2849 2850 !------------------------------------------------------------------------------!2851 ! Description:2852 ! ------------2853 !> Computation of sedimentation flux. Implementation according to Stevens2854 !> and Seifert (2008). Code is based on UCLA-LES.2855 !------------------------------------------------------------------------------!2856 SUBROUTINE sedimentation_rain2857 2858 IMPLICIT NONE2859 2860 INTEGER(iwp) :: i !< running index x direction2861 INTEGER(iwp) :: j !< running index y direction2862 INTEGER(iwp) :: k !< running index z direction2863 INTEGER(iwp) :: k_run !<2864 INTEGER(iwp) :: m !< running index surface elements2865 INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint2866 INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint2867 2868 REAL(wp) :: c_run !<2869 REAL(wp) :: d_max !<2870 REAL(wp) :: d_mean !<2871 REAL(wp) :: d_min !<2872 REAL(wp) :: dr !<2873 REAL(wp) :: flux !<2874 REAL(wp) :: flag !< flag to mask topography grid points2875 REAL(wp) :: lambda_r !<2876 REAL(wp) :: mu_r !<2877 REAL(wp) :: z_run !<2878 2879 REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !<2880 REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !<2881 REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !<2882 REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !<2883 REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nr !<2884 REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qr !<2885 REAL(wp), DIMENSION(nzb:nzt+1) :: w_nr !<2886 REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr !<2887 2888 CALL cpu_log( log_point_s(60), 'sed_rain', 'start' )2889 2890 !2891 !-- Compute velocities2892 DO i = nxlg, nxrg2893 DO j = nysg, nyng2894 DO k = nzb+1, nzt2895 !2896 !-- Predetermine flag to mask topography2897 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )2898 2899 IF ( qr(k,j,i) > eps_sb ) THEN2900 !2901 !-- Weight averaged diameter of rain drops:2902 dr = ( hyrho(k) * qr(k,j,i) / &2903 nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )2904 !2905 !-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005;2906 !-- Stevens and Seifert, 2008):2907 mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * &2908 ( dr - 1.4E-3_wp ) ) )2909 !2910 !-- Slope parameter of gamma distribution (Seifert, 2008):2911 lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * &2912 ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr2913 2914 w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &2915 a_term - b_term * ( 1.0_wp + &2916 c_term / &2917 lambda_r )**( -1.0_wp * &2918 ( mu_r + 1.0_wp ) ) &2919 ) &2920 ) * flag2921 2922 w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, &2923 a_term - b_term * ( 1.0_wp + &2924 c_term / &2925 lambda_r )**( -1.0_wp * &2926 ( mu_r + 4.0_wp ) ) &2927 ) &2928 ) * flag2929 ELSE2930 w_nr(k) = 0.0_wp2931 w_qr(k) = 0.0_wp2932 ENDIF2933 ENDDO2934 !2935 !-- Adjust boundary values using surface data type.2936 !-- Upward-facing2937 surf_s = bc_h(0)%start_index(j,i)2938 surf_e = bc_h(0)%end_index(j,i)2939 DO m = surf_s, surf_e2940 k = bc_h(0)%k(m)2941 w_nr(k-1) = w_nr(k)2942 w_qr(k-1) = w_qr(k)2943 ENDDO2944 !2945 !-- Downward-facing2946 surf_s = bc_h(1)%start_index(j,i)2947 surf_e = bc_h(1)%end_index(j,i)2948 DO m = surf_s, surf_e2949 k = bc_h(1)%k(m)2950 w_nr(k+1) = w_nr(k)2951 w_qr(k+1) = w_qr(k)2952 ENDDO2953 !2954 !-- Model top boundary value2955 w_nr(nzt+1) = 0.0_wp2956 w_qr(nzt+1) = 0.0_wp2957 !2958 !-- Compute Courant number2959 DO k = nzb+1, nzt2960 !2961 !-- Predetermine flag to mask topography2962 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )2963 2964 c_nr(k) = 0.25_wp * ( w_nr(k-1) + &2965 2.0_wp * w_nr(k) + w_nr(k+1) ) * &2966 dt_micro * ddzu(k) * flag2967 c_qr(k) = 0.25_wp * ( w_qr(k-1) + &2968 2.0_wp * w_qr(k) + w_qr(k+1) ) * &2969 dt_micro * ddzu(k) * flag2970 ENDDO2971 !2972 !-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977):2973 IF ( limiter_sedimentation ) THEN2974 2975 DO k = nzb+1, nzt2976 !2977 !-- Predetermine flag to mask topography2978 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )2979 2980 d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) )2981 d_min = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) )2982 d_max = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i)2983 2984 qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &2985 2.0_wp * d_max, &2986 ABS( d_mean ) ) &2987 * flag2988 2989 d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) )2990 d_min = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )2991 d_max = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)2992 2993 nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, &2994 2.0_wp * d_max, &2995 ABS( d_mean ) )2996 ENDDO2997 2998 ELSE2999 3000 nr_slope = 0.0_wp3001 qr_slope = 0.0_wp3002 3003 ENDIF3004 3005 sed_nr(nzt+1) = 0.0_wp3006 sed_qr(nzt+1) = 0.0_wp3007 !3008 !-- Compute sedimentation flux3009 DO k = nzt, nzb+1, -13010 !3011 !-- Predetermine flag to mask topography3012 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )3013 !3014 !-- Sum up all rain drop number densities which contribute to the flux3015 !-- through k-1/23016 flux = 0.0_wp3017 z_run = 0.0_wp ! height above z(k)3018 k_run = k3019 c_run = MIN( 1.0_wp, c_nr(k) )3020 DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )3021 flux = flux + hyrho(k_run) * &3022 ( nr(k_run,j,i) + nr_slope(k_run) * &3023 ( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run) &3024 * flag3025 z_run = z_run + dzu(k_run) * flag3026 k_run = k_run + 1 * flag3027 c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) ) &3028 * flag3029 ENDDO3030 !3031 !-- It is not allowed to sediment more rain drop number density than3032 !-- available3033 flux = MIN( flux, &3034 hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) * &3035 dt_micro &3036 )3037 3038 sed_nr(k) = flux / dt_micro * flag3039 nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) * &3040 ddzu(k+1) / hyrho(k) * dt_micro * flag3041 !3042 !-- Sum up all rain water content which contributes to the flux3043 !-- through k-1/23044 flux = 0.0_wp3045 z_run = 0.0_wp ! height above z(k)3046 k_run = k3047 c_run = MIN( 1.0_wp, c_qr(k) )3048 3049 DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt )3050 3051 flux = flux + hyrho(k_run) * ( qr(k_run,j,i) + &3052 qr_slope(k_run) * ( 1.0_wp - c_run ) * &3053 0.5_wp ) * c_run * dzu(k_run) * flag3054 z_run = z_run + dzu(k_run) * flag3055 k_run = k_run + 1 * flag3056 c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) ) &3057 * flag3058 3059 ENDDO3060 !3061 !-- It is not allowed to sediment more rain water content than3062 !-- available3063 flux = MIN( flux, &3064 hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) * &3065 dt_micro &3066 )3067 3068 sed_qr(k) = flux / dt_micro * flag3069 3070 qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * &3071 ddzu(k+1) / hyrho(k) * dt_micro * flag3072 q(k,j,i) = q(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * &3073 ddzu(k+1) / hyrho(k) * dt_micro * flag3074 pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) * &3075 ddzu(k+1) / hyrho(k) * lv_d_cp * &3076 d_exner(k) * dt_micro * flag3077 !3078 !-- Compute the rain rate3079 IF ( call_microphysics_at_all_substeps ) THEN3080 prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) &3081 * weight_substep(intermediate_timestep_count) &3082 * flag3083 ELSE3084 prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag3085 ENDIF3086 3087 ENDDO3088 ENDDO3089 ENDDO3090 3091 CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' )3092 3093 END SUBROUTINE sedimentation_rain3094 3095 3096 !------------------------------------------------------------------------------!3097 ! Description:3098 ! ------------3099 !> Computation of the precipitation amount due to gravitational settling of3100 !> rain and cloud (fog) droplets3101 !------------------------------------------------------------------------------!3102 SUBROUTINE calc_precipitation_amount3103 3104 IMPLICIT NONE3105 3106 INTEGER(iwp) :: i !< running index x direction3107 INTEGER(iwp) :: j !< running index y direction3108 INTEGER(iwp) :: k !< running index y direction3109 INTEGER(iwp) :: m !< running index surface elements3110 3111 IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.&3112 ( .NOT. call_microphysics_at_all_substeps .OR. &3113 intermediate_timestep_count == intermediate_timestep_count_max ) ) &3114 THEN3115 !3116 !-- Run over all upward-facing surface elements, i.e. non-natural,3117 !-- natural and urban3118 DO m = 1, bc_h(0)%ns3119 i = bc_h(0)%i(m)3120 j = bc_h(0)%j(m)3121 k = bc_h(0)%k(m)3122 precipitation_amount(j,i) = precipitation_amount(j,i) + &3123 prr(k,j,i) * hyrho(k) * dt_3d3124 ENDDO3125 3126 ENDIF3127 3128 END SUBROUTINE calc_precipitation_amount3129 3130 3131 !------------------------------------------------------------------------------!3132 ! Description:3133 ! ------------3134 !> Control of microphysics for grid points i,j3135 !------------------------------------------------------------------------------!3136 3137 SUBROUTINE bcm_actions_ij( i, j )3138 3139 IMPLICIT NONE3140 3141 INTEGER(iwp) :: i !<3142 INTEGER(iwp) :: j !<3143 3144 IF ( large_scale_forcing .AND. lsf_surf ) THEN3145 !3146 !-- Calculate vertical profile of the hydrostatic pressure (hyp)3147 hyp = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp)3148 d_exner = exner_function_invers(hyp)3149 exner = 1.0_wp / exner_function_invers(hyp)3150 hyrho = ideal_gas_law_rho_pt(hyp, pt_init)3151 !3152 !-- Compute reference density3153 rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp))3154 ENDIF3155 3156 !3157 !-- Compute length of time step3158 IF ( call_microphysics_at_all_substeps ) THEN3159 dt_micro = dt_3d * weight_pres(intermediate_timestep_count)3160 ELSE3161 dt_micro = dt_3d3162 ENDIF3163 !3164 !-- Reset precipitation rate3165 IF ( intermediate_timestep_count == 1 ) prr(:,j,i) = 0.0_wp3166 3167 !3168 !-- Compute cloud physics3169 IF( microphysics_kessler ) THEN3170 3171 CALL autoconversion_kessler( i,j )3172 IF ( cloud_water_sedimentation ) CALL sedimentation_cloud( i,j )3173 3174 ELSEIF ( microphysics_seifert ) THEN3175 3176 CALL adjust_cloud( i,j )3177 IF ( microphysics_morrison ) CALL activation( i,j )3178 IF ( microphysics_morrison ) CALL condensation( i,j )3179 CALL autoconversion( i,j )3180 CALL accretion( i,j )3181 CALL selfcollection_breakup( i,j )3182 CALL evaporation_rain( i,j )3183 CALL sedimentation_rain( i,j )3184 IF ( cloud_water_sedimentation ) CALL sedimentation_cloud( i,j )3185 3186 ENDIF3187 3188 CALL calc_precipitation_amount( i,j )3189 3190 END SUBROUTINE bcm_actions_ij3191 3192 !------------------------------------------------------------------------------!3193 ! Description:3194 ! ------------3195 !> Adjust number of raindrops to avoid nonlinear effects in3196 !> sedimentation and evaporation of rain drops due to too small or3197 !> too big weights of rain drops (Stevens and Seifert, 2008).3198 !> The same procedure is applied to cloud droplets if they are determined3199 !> prognostically. Call for grid point i,j3200 !------------------------------------------------------------------------------!3201 SUBROUTINE adjust_cloud_ij( i, j )3202 3203 IMPLICIT NONE3204 3205 INTEGER(iwp) :: i !<3206 INTEGER(iwp) :: j !<3207 INTEGER(iwp) :: k !<3208 3209 REAL(wp) :: flag !< flag to indicate first grid level above surface3210 3211 DO k = nzb+1, nzt3212 !3213 !-- Predetermine flag to mask topography3214 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )3215 3216 IF ( qr(k,j,i) <= eps_sb ) THEN3217 qr(k,j,i) = 0.0_wp3218 nr(k,j,i) = 0.0_wp3219 ELSE3220 !3221 !-- Adjust number of raindrops to avoid nonlinear effects in3222 !-- sedimentation and evaporation of rain drops due to too small or3223 !-- too big weights of rain drops (Stevens and Seifert, 2008).3224 IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) ) THEN3225 nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag3226 ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) ) THEN3227 nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag3228 ENDIF3229 3230 ENDIF3231 3232 IF ( microphysics_morrison ) THEN3233 IF ( qc(k,j,i) <= eps_sb ) THEN3234 qc(k,j,i) = 0.0_wp3235 nc(k,j,i) = 0.0_wp3236 ELSE3237 IF ( nc(k,j,i) * xcmin > qc(k,j,i) * hyrho(k) ) THEN3238 nc(k,j,i) = qc(k,j,i) * hyrho(k) / xamin * flag3239 ENDIF3240 ENDIF3241 ENDIF3242 3243 ENDDO3244 3245 END SUBROUTINE adjust_cloud_ij3246 3247 !------------------------------------------------------------------------------!3248 ! Description:3249 ! ------------3250 !> Calculate number of activated condensation nucleii after simple activation3251 !> scheme of Twomey, 1959.3252 !------------------------------------------------------------------------------!3253 SUBROUTINE activation_ij( i, j )3254 3255 IMPLICIT NONE3256 3257 INTEGER(iwp) :: i !<3258 INTEGER(iwp) :: j !<3259 INTEGER(iwp) :: k !<3260 3261 REAL(wp) :: activ !<3262 REAL(wp) :: afactor !<3263 REAL(wp) :: beta_act !<3264 REAL(wp) :: bfactor !<3265 REAL(wp) :: flag !< flag to indicate first grid level above surface3266 REAL(wp) :: k_act !<3267 REAL(wp) :: n_act !<3268 REAL(wp) :: n_ccn !<3269 REAL(wp) :: s_0 !<3270 REAL(wp) :: sat_max !<3271 REAL(wp) :: sigma !<3272 REAL(wp) :: sigma_act !<3273 3274 DO k = nzb+1, nzt3275 !3276 !-- Predetermine flag to mask topography3277 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )3278 !3279 !-- Call calculation of supersaturation3280 CALL supersaturation ( i, j, k )3281 !3282 !-- Prescribe parameters for activation3283 !-- (see: Bott + Trautmann, 2002, Atm. Res., 64)3284 k_act = 0.7_wp3285 activ = 0.0_wp3286 3287 IF ( sat > 0.0 .AND. .NOT. curvature_solution_effects_bulk ) THEN3288 !3289 !-- Compute the number of activated Aerosols3290 !-- (see: Twomey, 1959, Pure and applied Geophysics, 43)3291 n_act = na_init * sat**k_act3292 !3293 !-- Compute the number of cloud droplets3294 !-- (see: Morrison + Grabowski, 2007, JAS, 64)3295 ! activ = MAX( n_act - nc_d1(k), 0.0_wp) / dt_micro3296 3297 !3298 !-- Compute activation rate after Khairoutdinov and Kogan3299 !-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev., 128)3300 sat_max = 0.8_wp / 100.0_wp3301 activ = MAX( 0.0_wp, ( (na_init + nc(k,j,i) ) * MIN &3302 ( 1.0_wp, ( sat / sat_max )**k_act) - nc(k,j,i) ) ) / &3303 dt_micro3304 3305 nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro), na_init)3306 ELSEIF ( sat > 0.0 .AND. curvature_solution_effects_bulk ) THEN3307 !3308 !-- Curvature effect (afactor) with surface tension3309 !-- parameterization by Straka (2009)3310 sigma = 0.0761_wp - 0.000155_wp * ( t_l - 273.15_wp )3311 afactor = 2.0_wp * sigma / ( rho_l * r_v * t_l )3312 !3313 !-- Solute effect (bfactor)3314 bfactor = vanthoff * molecular_weight_of_water * &3315 rho_s / ( molecular_weight_of_solute * rho_l )3316 3317 !3318 !-- Prescribe power index that describes the soluble fraction3319 !-- of an aerosol particle (beta).3320 !-- (see: Morrison + Grabowski, 2007, JAS, 64)3321 beta_act = 0.5_wp3322 sigma_act = sigma_bulk**( 1.0_wp + beta_act )3323 !3324 !-- Calculate mean geometric supersaturation (s_0) with3325 !-- parameterization by Khvorostyanov and Curry (2006)3326 s_0 = dry_aerosol_radius **(- ( 1.0_wp + beta_act ) ) * &3327 ( 4.0_wp * afactor**3 / ( 27.0_wp * bfactor ) )**0.5_wp3328 3329 !3330 !-- Calculate number of activated CCN as a function of3331 !-- supersaturation and taking Koehler theory into account3332 !-- (see: Khvorostyanov + Curry, 2006, J. Geo. Res., 111)3333 n_ccn = ( na_init / 2.0_wp ) * ( 1.0_wp - ERF( &3334 LOG( s_0 / sat ) / ( SQRT(2.0_wp) * LOG(sigma_act) ) ) )3335 activ = MAX( ( n_ccn ) / dt_micro, 0.0_wp )3336 3337 nc(k,j,i) = MIN( (nc(k,j,i) + activ * dt_micro * flag), na_init)3338 ENDIF3339 3340 ENDDO3341 3342 END SUBROUTINE activation_ij3343 3344 !------------------------------------------------------------------------------!3345 ! Description:3346 ! ------------3347 !> Calculate condensation rate for cloud water content (after Khairoutdinov and3348 !> Kogan, 2000).3349 !------------------------------------------------------------------------------!3350 SUBROUTINE condensation_ij( i, j )3351 3352 IMPLICIT NONE3353 3354 INTEGER(iwp) :: i !<3355 INTEGER(iwp) :: j !<3356 INTEGER(iwp) :: k !<3357 3358 REAL(wp) :: cond !<3359 REAL(wp) :: cond_max !<3360 REAL(wp) :: dc !<3361 REAL(wp) :: evap !<3362 REAL(wp) :: flag !< flag to indicate first grid level above surface3363 REAL(wp) :: g_fac !<3364 REAL(wp) :: nc_0 !<3365 REAL(wp) :: temp !<3366 REAL(wp) :: xc !<3367 3368 3369 DO k = nzb+1, nzt3370 !3371 !-- Predetermine flag to mask topography3372 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )3373 !3374 !-- Call calculation of supersaturation3375 CALL supersaturation ( i, j, k )3376 !3377 !-- Actual temperature:3378 temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) )3379 3380 g_fac = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * &3381 l_v / ( thermal_conductivity_l * temp ) &3382 + r_v * temp / ( diff_coeff_l * e_s ) &3383 )3384 !3385 !-- Mean weight of cloud drops3386 IF ( nc(k,j,i) <= 0.0_wp) CYCLE3387 xc = MAX( (hyrho(k) * qc(k,j,i) / nc(k,j,i)), xcmin)3388 !3389 !-- Weight averaged diameter of cloud drops:3390 dc = ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )3391 !3392 !-- Integral diameter of cloud drops3393 nc_0 = nc(k,j,i) * dc3394 !3395 !-- Condensation needs only to be calculated in supersaturated regions3396 IF ( sat > 0.0_wp ) THEN3397 !3398 !-- Condensation rate of cloud water content3399 !-- after KK scheme.3400 !-- (see: Khairoutdinov + Kogan, 2000, Mon. Wea. Rev.,128)3401 cond = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)3402 cond_max = q(k,j,i) - q_s - qc(k,j,i) - qr(k,j,i)3403 cond = MIN( cond, cond_max / dt_micro )3404 3405 qc(k,j,i) = qc(k,j,i) + cond * dt_micro * flag3406 ELSEIF ( sat < 0.0_wp ) THEN3407 evap = 2.0_wp * pi * nc_0 * g_fac * sat / hyrho(k)3408 evap = MAX( evap, -qc(k,j,i) / dt_micro )3409 3410 qc(k,j,i) = qc(k,j,i) + evap * dt_micro * flag3411 ENDIF3412 ENDDO3413 3414 END SUBROUTINE condensation_ij3415 3416 3417 !------------------------------------------------------------------------------!3418 ! Description:3419 ! ------------3420 !> Autoconversion rate (Seifert and Beheng, 2006). Call for grid point i,j3421 !------------------------------------------------------------------------------!3422 SUBROUTINE autoconversion_ij( i, j )3423 3424 IMPLICIT NONE3425 3426 INTEGER(iwp) :: i !<3427 INTEGER(iwp) :: j !<3428 INTEGER(iwp) :: k !<3429 3430 REAL(wp) :: alpha_cc !<3431 REAL(wp) :: autocon !<3432 REAL(wp) :: dissipation !<3433 REAL(wp) :: flag !< flag to indicate first grid level above surface3434 REAL(wp) :: k_au !<3435 REAL(wp) :: l_mix !<3436 REAL(wp) :: nc_auto !<3437 REAL(wp) :: nu_c !<3438 REAL(wp) :: phi_au !<3439 REAL(wp) :: r_cc !<3440 REAL(wp) :: rc !<3441 REAL(wp) :: re_lambda !<3442 REAL(wp) :: sigma_cc !<3443 REAL(wp) :: tau_cloud !<3444 REAL(wp) :: xc !<3445 3446 DO k = nzb+1, nzt3447 !3448 !-- Predetermine flag to mask topography3449 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )3450 IF ( microphysics_morrison ) THEN3451 nc_auto = nc(k,j,i)3452 ELSE3453 nc_auto = nc_const3454 ENDIF3455 3456 IF ( qc(k,j,i) > eps_sb .AND. nc_auto > eps_mr ) THEN3457 3458 k_au = k_cc / ( 20.0_wp * x0 )3459 !3460 !-- Intern time scale of coagulation (Seifert and Beheng, 2006):3461 !-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ))3462 tau_cloud = MAX( 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) ), &3463 0.0_wp )3464 !3465 !-- Universal function for autoconversion process3466 !-- (Seifert and Beheng, 2006):3467 phi_au = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**33468 !3469 !-- Shape parameter of gamma distribution (Geoffroy et al., 2010):3470 !-- (Use constant nu_c = 1.0_wp instead?)3471 nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp )3472 !3473 !-- Mean weight of cloud droplets:3474 xc = hyrho(k) * qc(k,j,i) / nc_auto3475 !3476 !-- Parameterized turbulence effects on autoconversion (Seifert,3477 !-- Nuijens and Stevens, 2010)3478 IF ( collision_turbulence ) THEN3479 !3480 !-- Weight averaged radius of cloud droplets:3481 rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp )3482 3483 alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c )3484 r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c )3485 sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c )3486 !3487 !-- Mixing length (neglecting distance to ground and stratification)3488 l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp )3489 !3490 !-- Limit dissipation rate according to Seifert, Nuijens and3491 !-- Stevens (2010)3492 dissipation = MIN( 0.06_wp, diss(k,j,i) )3493 !3494 !-- Compute Taylor-microscale Reynolds number:3495 re_lambda = 6.0_wp / 11.0_wp * &3496 ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) * &3497 SQRT( 15.0_wp / kin_vis_air ) * &3498 dissipation**( 1.0_wp / 6.0_wp )3499 !3500 !-- The factor of 1.0E4 is needed to convert the dissipation rate3501 !-- from m2 s-3 to cm2 s-3.3502 k_au = k_au * ( 1.0_wp + &3503 dissipation * 1.0E4_wp * &3504 ( re_lambda * 1.0E-3_wp )**0.25_wp * &3505 ( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) / &3506 sigma_cc )**2 &3507 ) + beta_cc &3508 ) &3509 )3510 ENDIF3511 !3512 !-- Autoconversion rate (Seifert and Beheng, 2006):3513 autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / &3514 ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 * &3515 ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) * &3516 rho_surface3517 autocon = MIN( autocon, qc(k,j,i) / dt_micro )3518 3519 qr(k,j,i) = qr(k,j,i) + autocon * dt_micro * flag3520 qc(k,j,i) = qc(k,j,i) - autocon * dt_micro * flag3521 nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro * flag3522 IF ( microphysics_morrison ) THEN3523 nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), 2.0_wp * &3524 autocon / x0 * hyrho(k) * dt_micro * flag )3525 ENDIF3526 3527 ENDIF3528 3529 ENDDO3530 3531 END SUBROUTINE autoconversion_ij3532 3533 !------------------------------------------------------------------------------!3534 ! Description:3535 ! ------------3536 !> Autoconversion process (Kessler, 1969).3537 !------------------------------------------------------------------------------!3538 SUBROUTINE autoconversion_kessler_ij( i, j )3539 3540 3541 IMPLICIT NONE3542 3543 INTEGER(iwp) :: i !<3544 INTEGER(iwp) :: j !<3545 INTEGER(iwp) :: k !<3546 INTEGER(iwp) :: k_wall !< topography top index3547 3548 REAL(wp) :: dqdt_precip !<3549 REAL(wp) :: flag !< flag to indicate first grid level above surface3550 3551 !3552 !-- Determine vertical index of topography top3553 k_wall = get_topography_top_index_ji( j, i, 's' )3554 DO k = nzb+1, nzt3555 !3556 !-- Predetermine flag to mask topography3557 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )3558 3559 IF ( qc(k,j,i) > ql_crit ) THEN3560 dqdt_precip = prec_time_const * ( qc(k,j,i) - ql_crit )3561 ELSE3562 dqdt_precip = 0.0_wp3563 ENDIF3564 3565 qc(k,j,i) = qc(k,j,i) - dqdt_precip * dt_micro * flag3566 q(k,j,i) = q(k,j,i) - dqdt_precip * dt_micro * flag3567 pt(k,j,i) = pt(k,j,i) + dqdt_precip * dt_micro * lv_d_cp * d_exner(k) &3568 * flag3569 3570 !3571 !-- Compute the rain rate (stored on surface grid point)3572 prr(k_wall,j,i) = prr(k_wall,j,i) + dqdt_precip * dzw(k) * flag3573 3574 ENDDO3575 3576 END SUBROUTINE autoconversion_kessler_ij3577 3578 !------------------------------------------------------------------------------!3579 ! Description:3580 ! ------------3581 !> Accretion rate (Seifert and Beheng, 2006). Call for grid point i,j3582 !------------------------------------------------------------------------------!3583 SUBROUTINE accretion_ij( i, j )3584 3585 IMPLICIT NONE3586 3587 INTEGER(iwp) :: i !<3588 INTEGER(iwp) :: j !<3589 INTEGER(iwp) :: k !<3590 3591 REAL(wp) :: accr !<3592 REAL(wp) :: flag !< flag to indicate first grid level above surface3593 REAL(wp) :: k_cr !<3594 REAL(wp) :: nc_accr !<3595 REAL(wp) :: phi_ac !<3596 REAL(wp) :: tau_cloud !<3597 REAL(wp) :: xc !<3598 3599 3600 DO k = nzb+1, nzt3601 !3602 !-- Predetermine flag to mask topography3603 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )3604 IF ( microphysics_morrison ) THEN3605 nc_accr = nc(k,j,i)3606 ELSE3607 nc_accr = nc_const3608 ENDIF3609 3610 IF ( ( qc(k,j,i) > eps_sb ) .AND. ( qr(k,j,i) > eps_sb ) .AND. &3611 ( nc_accr > eps_mr ) ) THEN3612 !3613 !-- Intern time scale of coagulation (Seifert and Beheng, 2006):3614 tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )3615 !3616 !-- Universal function for accretion process3617 !-- (Seifert and Beheng, 2001):3618 phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**43619 3620 !3621 !-- Mean weight of cloud drops3622 xc = MAX( (hyrho(k) * qc(k,j,i) / nc_accr), xcmin)3623 !3624 !-- Parameterized turbulence effects on autoconversion (Seifert,3625 !-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to3626 !-- convert the dissipation rate (diss) from m2 s-3 to cm2 s-3.3627 IF ( collision_turbulence ) THEN3628 k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * &3629 MIN( 600.0_wp, &3630 diss(k,j,i) * 1.0E4_wp )**0.25_wp &3631 )3632 ELSE3633 k_cr = k_cr03634 ENDIF3635 !3636 !-- Accretion rate (Seifert and Beheng, 2006):3637 accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac * &3638 SQRT( rho_surface * hyrho(k) )3639 accr = MIN( accr, qc(k,j,i) / dt_micro )3640 3641 qr(k,j,i) = qr(k,j,i) + accr * dt_micro * flag3642 qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag3643 IF ( microphysics_morrison ) THEN3644 nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), accr / xc * &3645 hyrho(k) * dt_micro * flag &3646 )3647 ENDIF3648 3649 3650 ENDIF3651 3652 ENDDO3653 3654 END SUBROUTINE accretion_ij3655 3656 3657 !------------------------------------------------------------------------------!3658 ! Description:3659 ! ------------3660 !> Collisional breakup rate (Seifert, 2008). Call for grid point i,j3661 !------------------------------------------------------------------------------!3662 SUBROUTINE selfcollection_breakup_ij( i, j )3663 3664 IMPLICIT NONE3665 3666 INTEGER(iwp) :: i !<3667 INTEGER(iwp) :: j !<3668 INTEGER(iwp) :: k !<3669 3670 REAL(wp) :: breakup !<3671 REAL(wp) :: dr !<3672 REAL(wp) :: flag !< flag to indicate first grid level above surface3673 REAL(wp) :: phi_br !<3674 REAL(wp) :: selfcoll !<3675 3676 DO k = nzb+1, nzt3677 !3678 !-- Predetermine flag to mask topography3679 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )3680 3681 IF ( qr(k,j,i) > eps_sb ) THEN3682 !3683 !-- Selfcollection rate (Seifert and Beheng, 2001):3684 selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * SQRT( hyrho(k) * rho_surface )3685 !3686 !-- Weight averaged diameter of rain drops:3687 dr = ( hyrho(k) * qr(k,j,i) / nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp )3688 !3689 !-- Collisional breakup rate (Seifert, 2008):3690 IF ( dr >= 0.3E-3_wp ) THEN3691 phi_br = k_br * ( dr - 1.1E-3_wp )3692 breakup = selfcoll * ( phi_br + 1.0_wp )3693 ELSE3694 breakup = 0.0_wp3695 ENDIF3696 3697 selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro )3698 nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro * flag3699 3700 ENDIF3701 ENDDO3702 3703 END SUBROUTINE selfcollection_breakup_ij3704 3705 3706 !------------------------------------------------------------------------------!3707 ! Description:3708 ! ------------3709 3274 !> Evaporation of precipitable water. Condensation is neglected for 3710 3275 !> precipitable water. Call for grid point i,j … … 3739 3304 IF ( qr(k,j,i) > eps_sb ) THEN 3740 3305 ! 3741 !-- Call calculation of supersaturation 3306 !-- Call calculation of supersaturation 3742 3307 CALL supersaturation ( i, j, k ) 3743 3308 ! … … 3822 3387 ! ------------ 3823 3388 !> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR). 3389 !------------------------------------------------------------------------------! 3390 SUBROUTINE sedimentation_cloud 3391 3392 3393 IMPLICIT NONE 3394 3395 INTEGER(iwp) :: i !< 3396 INTEGER(iwp) :: j !< 3397 INTEGER(iwp) :: k !< 3398 3399 REAL(wp) :: flag !< flag to mask topography grid points 3400 REAL(wp) :: nc_sedi !< 3401 3402 REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !< 3403 REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nc !< 3404 3405 3406 CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' ) 3407 3408 sed_qc(nzt+1) = 0.0_wp 3409 sed_nc(nzt+1) = 0.0_wp 3410 3411 DO i = nxlg, nxrg 3412 DO j = nysg, nyng 3413 DO k = nzt, nzb+1, -1 3414 ! 3415 !-- Predetermine flag to mask topography 3416 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 3417 3418 IF ( microphysics_morrison ) THEN 3419 nc_sedi = nc(k,j,i) 3420 ELSE 3421 nc_sedi = nc_const 3422 ENDIF 3423 3424 ! 3425 !-- Sedimentation fluxes for number concentration are only calculated 3426 !-- for cloud_scheme = 'morrison' 3427 IF ( microphysics_morrison ) THEN 3428 IF ( qc(k,j,i) > eps_sb .AND. nc(k,j,i) > eps_mr ) THEN 3429 sed_nc(k) = sed_qc_const * & 3430 ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) * & 3431 ( nc(k,j,i) )**( 1.0_wp / 3.0_wp ) 3432 ELSE 3433 sed_nc(k) = 0.0_wp 3434 ENDIF 3435 3436 sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) * & 3437 nc(k,j,i) / dt_micro + sed_nc(k+1) & 3438 ) * flag 3439 3440 nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) * & 3441 ddzu(k+1) / hyrho(k) * dt_micro * flag 3442 ENDIF 3443 3444 IF ( qc(k,j,i) > eps_sb .AND. nc_sedi > eps_mr ) THEN 3445 sed_qc(k) = sed_qc_const * nc_sedi**( -2.0_wp / 3.0_wp ) * & 3446 ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp ) * & 3447 flag 3448 ELSE 3449 sed_qc(k) = 0.0_wp 3450 ENDIF 3451 3452 sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) / & 3453 dt_micro + sed_qc(k+1) & 3454 ) * flag 3455 3456 q(k,j,i) = q(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * & 3457 ddzu(k+1) / hyrho(k) * dt_micro * flag 3458 qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * & 3459 ddzu(k+1) / hyrho(k) * dt_micro * flag 3460 pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) * & 3461 ddzu(k+1) / hyrho(k) * lv_d_cp * & 3462 d_exner(k) * dt_micro * flag 3463 3464 ! 3465 !-- Compute the precipitation rate due to cloud (fog) droplets 3466 IF ( call_microphysics_at_all_substeps ) THEN 3467 prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) & 3468 * weight_substep(intermediate_timestep_count) & 3469 * flag 3470 ELSE 3471 prr(k,j,i) = prr(k,j,i) + sed_qc(k) / hyrho(k) * flag 3472 ENDIF 3473 3474 ENDDO 3475 ENDDO 3476 ENDDO 3477 3478 CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' ) 3479 3480 END SUBROUTINE sedimentation_cloud 3481 3482 3483 !------------------------------------------------------------------------------! 3484 ! Description: 3485 ! ------------ 3486 !> Sedimentation of cloud droplets (Ackermann et al., 2009, MWR). 3824 3487 !> Call for grid point i,j 3825 3488 !------------------------------------------------------------------------------! … … 3902 3565 3903 3566 END SUBROUTINE sedimentation_cloud_ij 3567 3568 3569 !------------------------------------------------------------------------------! 3570 ! Description: 3571 ! ------------ 3572 !> Computation of sedimentation flux. Implementation according to Stevens 3573 !> and Seifert (2008). Code is based on UCLA-LES. 3574 !------------------------------------------------------------------------------! 3575 SUBROUTINE sedimentation_rain 3576 3577 IMPLICIT NONE 3578 3579 INTEGER(iwp) :: i !< running index x direction 3580 INTEGER(iwp) :: j !< running index y direction 3581 INTEGER(iwp) :: k !< running index z direction 3582 INTEGER(iwp) :: k_run !< 3583 INTEGER(iwp) :: m !< running index surface elements 3584 INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint 3585 INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint 3586 3587 REAL(wp) :: c_run !< 3588 REAL(wp) :: d_max !< 3589 REAL(wp) :: d_mean !< 3590 REAL(wp) :: d_min !< 3591 REAL(wp) :: dr !< 3592 REAL(wp) :: flux !< 3593 REAL(wp) :: flag !< flag to mask topography grid points 3594 REAL(wp) :: lambda_r !< 3595 REAL(wp) :: mu_r !< 3596 REAL(wp) :: z_run !< 3597 3598 REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !< 3599 REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !< 3600 REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !< 3601 REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !< 3602 REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nr !< 3603 REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qr !< 3604 REAL(wp), DIMENSION(nzb:nzt+1) :: w_nr !< 3605 REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr !< 3606 3607 CALL cpu_log( log_point_s(60), 'sed_rain', 'start' ) 3608 3609 ! 3610 !-- Compute velocities 3611 DO i = nxlg, nxrg 3612 DO j = nysg, nyng 3613 DO k = nzb+1, nzt 3614 ! 3615 !-- Predetermine flag to mask topography 3616 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 3617 3618 IF ( qr(k,j,i) > eps_sb ) THEN 3619 ! 3620 !-- Weight averaged diameter of rain drops: 3621 dr = ( hyrho(k) * qr(k,j,i) / & 3622 nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp ) 3623 ! 3624 !-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005; 3625 !-- Stevens and Seifert, 2008): 3626 mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * & 3627 ( dr - 1.4E-3_wp ) ) ) 3628 ! 3629 !-- Slope parameter of gamma distribution (Seifert, 2008): 3630 lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * & 3631 ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr 3632 3633 w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, & 3634 a_term - b_term * ( 1.0_wp + & 3635 c_term / & 3636 lambda_r )**( -1.0_wp * & 3637 ( mu_r + 1.0_wp ) ) & 3638 ) & 3639 ) * flag 3640 3641 w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, & 3642 a_term - b_term * ( 1.0_wp + & 3643 c_term / & 3644 lambda_r )**( -1.0_wp * & 3645 ( mu_r + 4.0_wp ) ) & 3646 ) & 3647 ) * flag 3648 ELSE 3649 w_nr(k) = 0.0_wp 3650 w_qr(k) = 0.0_wp 3651 ENDIF 3652 ENDDO 3653 ! 3654 !-- Adjust boundary values using surface data type. 3655 !-- Upward-facing 3656 surf_s = bc_h(0)%start_index(j,i) 3657 surf_e = bc_h(0)%end_index(j,i) 3658 DO m = surf_s, surf_e 3659 k = bc_h(0)%k(m) 3660 w_nr(k-1) = w_nr(k) 3661 w_qr(k-1) = w_qr(k) 3662 ENDDO 3663 ! 3664 !-- Downward-facing 3665 surf_s = bc_h(1)%start_index(j,i) 3666 surf_e = bc_h(1)%end_index(j,i) 3667 DO m = surf_s, surf_e 3668 k = bc_h(1)%k(m) 3669 w_nr(k+1) = w_nr(k) 3670 w_qr(k+1) = w_qr(k) 3671 ENDDO 3672 ! 3673 !-- Model top boundary value 3674 w_nr(nzt+1) = 0.0_wp 3675 w_qr(nzt+1) = 0.0_wp 3676 ! 3677 !-- Compute Courant number 3678 DO k = nzb+1, nzt 3679 ! 3680 !-- Predetermine flag to mask topography 3681 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 3682 3683 c_nr(k) = 0.25_wp * ( w_nr(k-1) + & 3684 2.0_wp * w_nr(k) + w_nr(k+1) ) * & 3685 dt_micro * ddzu(k) * flag 3686 c_qr(k) = 0.25_wp * ( w_qr(k-1) + & 3687 2.0_wp * w_qr(k) + w_qr(k+1) ) * & 3688 dt_micro * ddzu(k) * flag 3689 ENDDO 3690 ! 3691 !-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977): 3692 IF ( limiter_sedimentation ) THEN 3693 3694 DO k = nzb+1, nzt 3695 ! 3696 !-- Predetermine flag to mask topography 3697 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 3698 3699 d_mean = 0.5_wp * ( qr(k+1,j,i) - qr(k-1,j,i) ) 3700 d_min = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) 3701 d_max = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i) 3702 3703 qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, & 3704 2.0_wp * d_max, & 3705 ABS( d_mean ) ) & 3706 * flag 3707 3708 d_mean = 0.5_wp * ( nr(k+1,j,i) - nr(k-1,j,i) ) 3709 d_min = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) 3710 d_max = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i) 3711 3712 nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, & 3713 2.0_wp * d_max, & 3714 ABS( d_mean ) ) 3715 ENDDO 3716 3717 ELSE 3718 3719 nr_slope = 0.0_wp 3720 qr_slope = 0.0_wp 3721 3722 ENDIF 3723 3724 sed_nr(nzt+1) = 0.0_wp 3725 sed_qr(nzt+1) = 0.0_wp 3726 ! 3727 !-- Compute sedimentation flux 3728 DO k = nzt, nzb+1, -1 3729 ! 3730 !-- Predetermine flag to mask topography 3731 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) ) 3732 ! 3733 !-- Sum up all rain drop number densities which contribute to the flux 3734 !-- through k-1/2 3735 flux = 0.0_wp 3736 z_run = 0.0_wp ! height above z(k) 3737 k_run = k 3738 c_run = MIN( 1.0_wp, c_nr(k) ) 3739 DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt ) 3740 flux = flux + hyrho(k_run) * & 3741 ( nr(k_run,j,i) + nr_slope(k_run) * & 3742 ( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run) & 3743 * flag 3744 z_run = z_run + dzu(k_run) * flag 3745 k_run = k_run + 1 * flag 3746 c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) ) & 3747 * flag 3748 ENDDO 3749 ! 3750 !-- It is not allowed to sediment more rain drop number density than 3751 !-- available 3752 flux = MIN( flux, & 3753 hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) * & 3754 dt_micro & 3755 ) 3756 3757 sed_nr(k) = flux / dt_micro * flag 3758 nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) * & 3759 ddzu(k+1) / hyrho(k) * dt_micro * flag 3760 ! 3761 !-- Sum up all rain water content which contributes to the flux 3762 !-- through k-1/2 3763 flux = 0.0_wp 3764 z_run = 0.0_wp ! height above z(k) 3765 k_run = k 3766 c_run = MIN( 1.0_wp, c_qr(k) ) 3767 3768 DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt ) 3769 3770 flux = flux + hyrho(k_run) * ( qr(k_run,j,i) + & 3771 qr_slope(k_run) * ( 1.0_wp - c_run ) * & 3772 0.5_wp ) * c_run * dzu(k_run) * flag 3773 z_run = z_run + dzu(k_run) * flag 3774 k_run = k_run + 1 * flag 3775 c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) ) & 3776 * flag 3777 3778 ENDDO 3779 ! 3780 !-- It is not allowed to sediment more rain water content than 3781 !-- available 3782 flux = MIN( flux, & 3783 hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) * & 3784 dt_micro & 3785 ) 3786 3787 sed_qr(k) = flux / dt_micro * flag 3788 3789 qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * & 3790 ddzu(k+1) / hyrho(k) * dt_micro * flag 3791 q(k,j,i) = q(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * & 3792 ddzu(k+1) / hyrho(k) * dt_micro * flag 3793 pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) * & 3794 ddzu(k+1) / hyrho(k) * lv_d_cp * & 3795 d_exner(k) * dt_micro * flag 3796 ! 3797 !-- Compute the rain rate 3798 IF ( call_microphysics_at_all_substeps ) THEN 3799 prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) & 3800 * weight_substep(intermediate_timestep_count) & 3801 * flag 3802 ELSE 3803 prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * flag 3804 ENDIF 3805 3806 ENDDO 3807 ENDDO 3808 ENDDO 3809 3810 CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' ) 3811 3812 END SUBROUTINE sedimentation_rain 3904 3813 3905 3814 … … 4123 4032 4124 4033 END SUBROUTINE sedimentation_rain_ij 4034 4035 4036 !------------------------------------------------------------------------------! 4037 ! Description: 4038 ! ------------ 4039 !> Computation of the precipitation amount due to gravitational settling of 4040 !> rain and cloud (fog) droplets 4041 !------------------------------------------------------------------------------! 4042 SUBROUTINE calc_precipitation_amount 4043 4044 IMPLICIT NONE 4045 4046 INTEGER(iwp) :: i !< running index x direction 4047 INTEGER(iwp) :: j !< running index y direction 4048 INTEGER(iwp) :: k !< running index y direction 4049 INTEGER(iwp) :: m !< running index surface elements 4050 4051 IF ( ( dt_do2d_xy - time_do2d_xy ) < precipitation_amount_interval .AND.& 4052 ( .NOT. call_microphysics_at_all_substeps .OR. & 4053 intermediate_timestep_count == intermediate_timestep_count_max ) ) & 4054 THEN 4055 ! 4056 !-- Run over all upward-facing surface elements, i.e. non-natural, 4057 !-- natural and urban 4058 DO m = 1, bc_h(0)%ns 4059 i = bc_h(0)%i(m) 4060 j = bc_h(0)%j(m) 4061 k = bc_h(0)%k(m) 4062 precipitation_amount(j,i) = precipitation_amount(j,i) + & 4063 prr(k,j,i) * hyrho(k) * dt_3d 4064 ENDDO 4065 4066 ENDIF 4067 4068 END SUBROUTINE calc_precipitation_amount 4125 4069 4126 4070
Note: See TracChangeset
for help on using the changeset viewer.
