Changeset 4502 for palm/trunk/SOURCE/bulk_cloud_model_mod.f90
- Timestamp:
- Apr 17, 2020 4:14:16 PM (4 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/bulk_cloud_model_mod.f90 (modified) (100 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/bulk_cloud_model_mod.f90
r4495 r4502 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Implementation of ice microphysics 28 ! 29 ! 4495 2020-04-13 20:11:20Z raasch 27 30 ! restart data handling with MPI-IO added 28 31 ! … … 116 119 precipitation_amount, prr, pt, d_exner, pt_init, q, ql, ql_1, & 117 120 qc, qc_1, qc_2, qc_3, qc_p, qr, qr_1, qr_2, qr_3, qr_p, & 118 exner, zu, tnc_m, tnr_m, tqc_m, tqr_m, tend, rdf_sc, & 119 flux_l_qc, flux_l_qr, flux_l_nc, flux_l_nr, & 120 flux_s_qc, flux_s_qr, flux_s_nc, flux_s_nr, & 121 diss_l_qc, diss_l_qr, diss_l_nc, diss_l_nr, & 122 diss_s_qc, diss_s_qr, diss_s_nc, diss_s_nr 121 exner, zu, tnc_m, tnr_m, tqc_m, tqr_m, tend, rdf_sc, & 122 flux_l_qc, flux_l_qr, flux_l_nc, flux_l_nr, & 123 flux_s_qc, flux_s_qr, flux_s_nc, flux_s_nr, & 124 diss_l_qc, diss_l_qr, diss_l_nc, diss_l_nr, & 125 diss_s_qc, diss_s_qr, diss_s_nc, diss_s_nr, & 126 ni, ni_1, ni_2, ni_3, ni_p, tni_m, & 127 qi, qi_1, qi_2, qi_3, qi_p, tqi_m, & 128 flux_l_qi, flux_l_ni, flux_s_qi, flux_s_ni, & 129 diss_l_qi, diss_l_ni, diss_s_qi, diss_s_ni 130 123 131 124 132 USE averaging, & 125 ONLY: nc_av, nr_av, prr_av, qc_av, ql_av, qr_av 133 ONLY: nc_av, nr_av, prr_av, qc_av, ql_av, qr_av, ni_av, qi_av 126 134 127 135 USE basic_constants_and_equations_mod, & 128 ONLY: c_p, g, lv_d_cp, lv_d_rd, l_v, magnus, molecular_weight_of_solute,& 136 ONLY: c_p, g, lv_d_cp, lv_d_rd, l_v, magnus, magnus_ice, & 137 molecular_weight_of_solute, & 129 138 molecular_weight_of_water, pi, rho_l, rho_s, r_d, r_v, vanthoff,& 130 139 exner_function, exner_function_invers, ideal_gas_law_rho, & 131 ideal_gas_law_rho_pt, barometric_formula, rd_d_rv 140 ideal_gas_law_rho_pt, barometric_formula, rd_d_rv, l_s, & 141 ls_d_cp 132 142 133 143 USE control_parameters, & … … 143 153 dt_3d, dt_do2d_xy, intermediate_timestep_count, & 144 154 intermediate_timestep_count_max, large_scale_forcing, & 145 lsf_surf, pt_surface, restart_data_format_output, rho_surface, surface_pressure, & 155 lsf_surf, pt_surface, restart_data_format_output, rho_surface, & 156 surface_pressure, & 146 157 time_do2d_xy, message_string, initializing_actions, & 147 ws_scheme_sca, scalar_advec, timestep_scheme, tsc, loop_optimization 158 ws_scheme_sca, scalar_advec, timestep_scheme, tsc, & 159 loop_optimization, simulated_time 148 160 149 161 USE cpulog, & … … 167 179 ONLY: threads_per_task 168 180 169 USE restart_data_mpi_io_mod, 181 USE restart_data_mpi_io_mod, & 170 182 ONLY: rrd_mpi_io, wrd_mpi_io 171 183 172 184 USE statistics, & 173 185 ONLY: weight_pres, weight_substep, sums_wsncs_ws_l, sums_wsnrs_ws_l, & 174 sums_wsqcs_ws_l, sums_wsqrs_ws_l 186 sums_wsqcs_ws_l, sums_wsqrs_ws_l, & 187 sums_wsqis_ws_l, sums_wsnis_ws_l 175 188 176 189 USE surface_mod, & 177 190 ONLY : bc_h, & 178 191 surf_bulk_cloud_model, & 179 surf_microphysics_morrison, surf_microphysics_seifert, & 180 surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, surf_usm_v 192 surf_microphysics_morrison, surf_microphysics_seifert, & 193 surf_microphysics_ice_extension, & 194 surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, & 195 surf_usm_v 181 196 182 197 IMPLICIT NONE … … 193 208 LOGICAL :: cloud_water_sedimentation = .FALSE. !< cloud water sedimentation 194 209 LOGICAL :: curvature_solution_effects_bulk = .FALSE. !< flag for considering koehler theory 210 LOGICAL :: ice_crystal_sedimentation = .FALSE. !< flag for ice crystal sedimentation 195 211 LOGICAL :: limiter_sedimentation = .TRUE. !< sedimentation limiter 196 212 LOGICAL :: collision_turbulence = .FALSE. !< turbulence effects … … 198 214 199 215 LOGICAL :: call_microphysics_at_all_substeps = .FALSE. !< namelist parameter 216 LOGICAL :: microphysics_ice_extension = .FALSE. !< use ice microphysics scheme 200 217 LOGICAL :: microphysics_sat_adjust = .FALSE. !< use saturation adjust bulk scheme? 201 218 LOGICAL :: microphysics_kessler = .FALSE. !< use kessler bulk scheme? 202 219 LOGICAL :: microphysics_morrison = .FALSE. !< use 2-moment Morrison (add. prog. eq. for nc and qc) 203 220 LOGICAL :: microphysics_seifert = .FALSE. !< use 2-moment Seifert and Beheng scheme 204 LOGICAL :: microphysics_morrison_no_rain = .FALSE. !< use 2-moment Morrison 221 LOGICAL :: microphysics_morrison_no_rain = .FALSE. !< use 2-moment Morrison 205 222 LOGICAL :: precipitation = .FALSE. !< namelist parameter 206 223 … … 240 257 REAL(wp) :: w_precipitation = 9.65_wp !< maximum terminal velocity (m s-1) 241 258 REAL(wp) :: x0 = 2.6E-10_wp !< separating drop mass (kg) 242 ! REAL(wp) :: xamin = 5.24E-19_wp !< average aerosol mass (kg) (~ 0.05µm)259 REAL(wp) :: ximin = 4.42E-14_wp !< minimum mass of ice crystal (kg) (D~10µm) 243 260 REAL(wp) :: xcmin = 4.18E-15_wp !< minimum cloud drop size (kg) (~ 1µm) 244 261 REAL(wp) :: xrmin = 2.6E-10_wp !< minimum rain drop size (kg) … … 249 266 REAL(wp) :: dry_aerosol_radius = 0.05E-6_wp !< dry aerosol radius 250 267 REAL(wp) :: dt_micro !< microphysics time step 268 REAL(wp) :: in_init = 1000.0_wp !< initial number of potential ice nucleii 251 269 REAL(wp) :: sigma_bulk = 2.0_wp !< width of aerosol spectrum 252 270 REAL(wp) :: na_init = 100.0E6_wp !< Total particle/aerosol concentration (cm-3) … … 254 272 REAL(wp) :: dt_precipitation = 100.0_wp !< timestep precipitation (s) 255 273 REAL(wp) :: sed_qc_const !< const. for sedimentation of cloud water 256 REAL(wp) :: pirho_l !< pi * rho_l / 6.0; 274 REAL(wp) :: pirho_l !< pi * rho_l / 6.0 275 REAL(wp) :: start_ice_microphysics = 0.0_wp !< time from which on ice microhysics are executed 257 276 258 277 REAL(wp) :: e_s !< saturation water vapor pressure 278 REAL(wp) :: e_si !< saturation water vapor pressure over ice 259 279 REAL(wp) :: q_s !< saturation mixing ratio 280 REAL(wp) :: q_si !< saturation mixing ratio over ice 260 281 REAL(wp) :: sat !< supersaturation 261 REAL(wp) :: t_l !< actual temperature 282 REAL(wp) :: sat_ice !< supersaturation over ice 283 REAL(wp) :: t_l !< liquid-(ice) water temperature 262 284 263 285 SAVE … … 294 316 dt_precipitation, & 295 317 microphysics_morrison, & 296 microphysics_morrison_no_rain, & 318 microphysics_morrison_no_rain, & 297 319 microphysics_sat_adjust, & 298 320 microphysics_seifert, & 321 microphysics_ice_extension, & 299 322 na_init, & 300 323 nc_const, & 301 324 precipitation, & 302 sigma_gc 303 325 sigma_gc, & 326 start_ice_microphysics, & 327 ice_crystal_sedimentation, & 328 in_init 304 329 305 330 INTERFACE bcm_parin … … 420 445 nc_const, & 421 446 sigma_bulk, & 422 ventilation_effect 447 ventilation_effect, & 448 ice_crystal_sedimentation, & 449 microphysics_ice_extension, & 450 start_ice_microphysics, & 451 in_init 423 452 424 453 line = ' ' … … 517 546 microphysics_kessler = .FALSE. 518 547 microphysics_morrison = .TRUE. 519 microphysics_morrison_no_rain = .TRUE. 520 precipitation = .FALSE. 548 microphysics_morrison_no_rain = .TRUE. 549 precipitation = .FALSE. 521 550 ELSE 522 551 message_string = 'unknown cloud microphysics scheme cloud_scheme ="' // & … … 546 575 surf_microphysics_morrison = microphysics_morrison 547 576 surf_microphysics_seifert = microphysics_seifert 548 577 surf_microphysics_ice_extension = microphysics_ice_extension 549 578 ! 550 579 !-- Check aerosol … … 592 621 unit = '1/m3' 593 622 623 CASE ( 'ni' ) 624 IF ( .NOT. microphysics_ice_extension ) THEN 625 message_string = 'output of "' // TRIM( var ) // '" ' // & 626 'requires ' // & 627 'microphysics_ice_extension = ".TRUE."' 628 CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 ) 629 ENDIF 630 unit = '1/m3' 631 594 632 CASE ( 'nr' ) 595 633 IF ( .NOT. microphysics_seifert ) THEN … … 611 649 612 650 CASE ( 'qc' ) 651 unit = 'kg/kg' 652 653 CASE ( 'qi' ) 654 IF ( .NOT. microphysics_ice_extension ) THEN 655 message_string = 'output of "' // TRIM( var ) // '" ' // & 656 'requires ' // & 657 'microphysics_ice_extension = ".TRUE."' 658 CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 ) 659 ENDIF 613 660 unit = 'kg/kg' 614 661 … … 698 745 ENDIF 699 746 pr_index = 123 747 dopr_index(var_count) = pr_index 748 dopr_unit = '1/m3' 749 unit = dopr_unit 750 hom(:,2,pr_index,:) = SPREAD( zu, 2, statistic_regions+1 ) 751 752 CASE ( 'ni' ) 753 IF ( .NOT. microphysics_ice_extension ) THEN 754 message_string = 'data_output_pr = ' // & 755 TRIM( data_output_pr(var_count) ) // & 756 ' is not implemented for' // & 757 ' microphysics_ice_extension = ".F."' 758 CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 ) 759 ENDIF 760 pr_index = 124 700 761 dopr_index(var_count) = pr_index 701 762 dopr_unit = '1/m3' … … 730 791 unit = dopr_unit 731 792 hom(:,2,pr_index,:) = SPREAD( zu, 2, statistic_regions+1 ) 793 732 794 CASE ( 'qc' ) 733 795 pr_index = 75 796 dopr_index(var_count) = pr_index 797 dopr_unit = 'kg/kg' 798 unit = dopr_unit 799 hom(:,2,pr_index,:) = SPREAD( zu, 2, statistic_regions+1 ) 800 801 CASE ( 'qi' ) 802 IF ( .NOT. microphysics_ice_extension ) THEN 803 message_string = 'data_output_pr = ' // & 804 TRIM( data_output_pr(var_count) ) // & 805 ' is not implemented for' // & 806 ' microphysics_ice_extension = ".F."' 807 CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 ) 808 ENDIF 809 pr_index = 125 734 810 dopr_index(var_count) = pr_index 735 811 dopr_unit = 'kg/kg' … … 792 868 ! 793 869 !-- 3D-cloud drop water content, cloud drop concentration arrays 794 ALLOCATE( nc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &795 nc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &796 nc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &797 qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &798 qc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &870 ALLOCATE( nc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 871 nc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 872 nc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 873 qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 874 qc_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 799 875 qc_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 800 876 ENDIF … … 803 879 ! 804 880 !-- 3D-rain water content, rain drop concentration arrays 805 ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &806 nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &807 nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &808 qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &809 qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &881 ALLOCATE( nr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 882 nr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 883 nr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 884 qr_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 885 qr_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 810 886 qr_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 811 887 ENDIF 812 888 889 IF ( microphysics_ice_extension ) THEN 890 ! 891 !-- 3D-cloud drop water content, cloud drop concentration arrays 892 ALLOCATE( ni_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 893 ni_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 894 ni_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 895 qi_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 896 qi_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 897 qi_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 898 ENDIF 899 813 900 IF ( ws_scheme_sca ) THEN 814 815 901 IF ( microphysics_morrison ) THEN 816 902 ALLOCATE( sums_wsqcs_ws_l(nzb:nzt+1,0:threads_per_task-1) ) … … 819 905 sums_wsncs_ws_l = 0.0_wp 820 906 ENDIF 821 822 907 IF ( microphysics_seifert ) THEN 823 908 ALLOCATE( sums_wsqrs_ws_l(nzb:nzt+1,0:threads_per_task-1) ) … … 826 911 sums_wsnrs_ws_l = 0.0_wp 827 912 ENDIF 828 913 IF ( microphysics_ice_extension ) THEN 914 ALLOCATE( sums_wsqis_ws_l(nzb:nzt+1,0:threads_per_task-1) ) 915 ALLOCATE( sums_wsnis_ws_l(nzb:nzt+1,0:threads_per_task-1) ) 916 sums_wsqis_ws_l = 0.0_wp 917 sums_wsnis_ws_l = 0.0_wp 918 ENDIF 829 919 ENDIF 830 920 … … 835 925 !-- advection routines. 836 926 IF ( loop_optimization /= 'vector' ) THEN 837 838 927 IF ( ws_scheme_sca ) THEN 839 840 928 IF ( microphysics_morrison ) THEN 841 929 ALLOCATE( flux_s_qc(nzb+1:nzt,0:threads_per_task-1), & … … 848 936 diss_l_nc(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) 849 937 ENDIF 850 851 938 IF ( microphysics_seifert ) THEN 852 939 ALLOCATE( flux_s_qr(nzb+1:nzt,0:threads_per_task-1), & … … 859 946 diss_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) 860 947 ENDIF 861 862 ENDIF 863 948 IF ( microphysics_ice_extension ) THEN 949 ALLOCATE( flux_s_qi(nzb+1:nzt,0:threads_per_task-1), & 950 diss_s_qi(nzb+1:nzt,0:threads_per_task-1), & 951 flux_s_ni(nzb+1:nzt,0:threads_per_task-1), & 952 diss_s_ni(nzb+1:nzt,0:threads_per_task-1) ) 953 ALLOCATE( flux_l_qi(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 954 diss_l_qi(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 955 flux_l_ni(nzb+1:nzt,nys:nyn,0:threads_per_task-1), & 956 diss_l_ni(nzb+1:nzt,nys:nyn,0:threads_per_task-1) ) 957 ENDIF 958 ENDIF 864 959 ENDIF 865 960 … … 878 973 nr => nr_1; nr_p => nr_2; tnr_m => nr_3 879 974 ENDIF 975 IF ( microphysics_ice_extension ) THEN 976 qi => qi_1; qi_p => qi_2; tqi_m => qi_3 977 ni => ni_1; ni_p => ni_2; tni_m => ni_3 978 ENDIF 979 880 980 881 981 … … 919 1019 ENDIF 920 1020 ! 1021 !-- Initialize the remaining quantities 1022 IF ( microphysics_ice_extension ) THEN 1023 DO i = nxlg, nxrg 1024 DO j = nysg, nyng 1025 qi(:,j,i) = 0.0_wp 1026 ni(:,j,i) = 0.0_wp 1027 ENDDO 1028 ENDDO 1029 ENDIF 1030 ! 921 1031 !-- Liquid water content and precipitation amount 922 1032 !-- are zero at beginning of the simulation … … 938 1048 qr_p = qr 939 1049 nr_p = nr 1050 ENDIF 1051 IF ( microphysics_ice_extension ) THEN 1052 tqi_m = 0.0_wp 1053 tni_m = 0.0_wp 1054 qi_p = qi 1055 ni_p = ni 940 1056 ENDIF 941 1057 ENDIF ! Only if not read_restart_data … … 1080 1196 sums_wsnrs_ws_l = 0.0_wp 1081 1197 ENDIF 1198 IF ( microphysics_ice_extension ) THEN 1199 sums_wsqis_ws_l = 0.0_wp 1200 sums_wsnis_ws_l = 0.0_wp 1201 ENDIF 1082 1202 1083 1203 ENDIF … … 1096 1216 !> Control of microphysics for grid points i,j 1097 1217 !------------------------------------------------------------------------------! 1098 1099 1218 SUBROUTINE bcm_actions_ij( i, j, location ) 1100 1219 … … 1121 1240 sums_wsnrs_ws_l = 0.0_wp 1122 1241 ENDIF 1242 IF ( microphysics_ice_extension ) THEN 1243 sums_wsqis_ws_l = 0.0_wp 1244 sums_wsnis_ws_l = 0.0_wp 1245 ENDIF 1246 1123 1247 1124 1248 ENDIF … … 1143 1267 CALL cpu_log( log_point(51), 'microphysics', 'start' ) 1144 1268 1145 IF ( .NOT. microphysics_sat_adjust .AND. &1146 ( intermediate_timestep_count == 1 .OR. 1147 call_microphysics_at_all_substeps ) ) 1269 IF ( .NOT. microphysics_sat_adjust .AND. & 1270 ( intermediate_timestep_count == 1 .OR. & 1271 call_microphysics_at_all_substeps ) ) & 1148 1272 THEN 1149 1273 … … 1151 1275 ! 1152 1276 !-- Calculate vertical profile of the hydrostatic pressure (hyp) 1153 hyp = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp) 1277 hyp = barometric_formula(zu, pt_surface * & 1278 exner_function(surface_pressure * 100.0_wp), & 1279 surface_pressure * 100.0_wp) 1154 1280 d_exner = exner_function_invers(hyp) 1155 1281 exner = 1.0_wp / exner_function_invers(hyp) 1156 hyrho = ideal_gas_law_rho_pt(hyp, pt _init)1282 hyrho = ideal_gas_law_rho_pt(hyp, pt(:, nys, nxl) ) 1157 1283 ! 1158 1284 !-- Compute reference density 1159 rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp)) 1285 rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, & 1286 pt_surface * & 1287 exner_function(surface_pressure * 100.0_wp)) 1160 1288 ENDIF 1161 1289 … … 1178 1306 CALL autoconversion_kessler 1179 1307 IF ( cloud_water_sedimentation ) CALL sedimentation_cloud 1180 1308 1181 1309 ! 1182 1310 !-- Here the seifert beheng scheme is used. Cloud concentration is assumed to … … 1184 1312 ELSEIF ( microphysics_seifert .AND. .NOT. microphysics_morrison ) THEN 1185 1313 CALL adjust_cloud 1314 IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics ) & 1315 CALL ice_nucleation 1316 IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics ) & 1317 CALL ice_deposition 1186 1318 CALL autoconversion 1187 1319 CALL accretion … … 1190 1322 CALL sedimentation_rain 1191 1323 IF ( cloud_water_sedimentation ) CALL sedimentation_cloud 1192 1193 ! 1194 !-- Here the morrison scheme is used. No rain processes are considered and qr and nr 1324 IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics ) & 1325 CALL adjust_ice 1326 IF ( ice_crystal_sedimentation .AND. microphysics_ice_extension & 1327 .AND. simulated_time > start_ice_microphysics ) CALL sedimentation_ice 1328 1329 ! 1330 !-- Here the morrison scheme is used. No rain processes are considered and qr and nr 1195 1331 !-- are not allocated 1196 1332 ELSEIF ( microphysics_morrison_no_rain .AND. .NOT. microphysics_seifert ) THEN 1197 1333 CALL activation 1198 1334 CALL condensation 1199 IF ( cloud_water_sedimentation ) CALL sedimentation_cloud 1200 1201 ! 1202 !-- Here the full morrison scheme is used and all processes of Seifert and Beheng are 1335 IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics ) & 1336 CALL adjust_ice 1337 IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics ) & 1338 CALL ice_nucleation 1339 IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics ) & 1340 CALL ice_deposition 1341 IF ( cloud_water_sedimentation ) CALL sedimentation_cloud 1342 1343 ! 1344 !-- Here the full morrison scheme is used and all processes of Seifert and Beheng are 1203 1345 !-- included 1204 1346 ELSEIF ( microphysics_morrison .AND. microphysics_seifert ) THEN … … 1206 1348 CALL activation 1207 1349 CALL condensation 1350 IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics ) & 1351 CALL adjust_ice 1352 IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics ) & 1353 CALL ice_nucleation 1354 IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics ) & 1355 CALL ice_deposition 1208 1356 CALL autoconversion 1209 1357 CALL accretion … … 1211 1359 CALL evaporation_rain 1212 1360 CALL sedimentation_rain 1213 IF ( cloud_water_sedimentation ) CALL sedimentation_cloud 1361 IF ( cloud_water_sedimentation ) CALL sedimentation_cloud 1214 1362 1215 1363 ENDIF … … 1229 1377 !> Control of microphysics for grid points i,j 1230 1378 !------------------------------------------------------------------------------! 1231 1232 1379 SUBROUTINE bcm_non_advective_processes_ij( i, j ) 1233 1380 … … 1236 1383 INTEGER(iwp) :: j !< 1237 1384 1238 IF ( .NOT. microphysics_sat_adjust .AND. &1239 ( intermediate_timestep_count == 1 .OR. 1240 call_microphysics_at_all_substeps ) ) 1385 IF ( .NOT. microphysics_sat_adjust .AND. & 1386 ( intermediate_timestep_count == 1 .OR. & 1387 call_microphysics_at_all_substeps ) ) & 1241 1388 THEN 1242 1389 … … 1244 1391 ! 1245 1392 !-- Calculate vertical profile of the hydrostatic pressure (hyp) 1246 hyp = barometric_formula(zu, pt_surface * exner_function(surface_pressure * 100.0_wp), surface_pressure * 100.0_wp) 1393 hyp = barometric_formula(zu, pt_surface * & 1394 exner_function(surface_pressure * 100.0_wp), & 1395 surface_pressure * 100.0_wp) 1247 1396 d_exner = exner_function_invers(hyp) 1248 1397 exner = 1.0_wp / exner_function_invers(hyp) 1249 hyrho = ideal_gas_law_rho_pt(hyp, pt _init)1398 hyrho = ideal_gas_law_rho_pt(hyp, pt(:, nys, nxl) ) 1250 1399 ! 1251 1400 !-- Compute reference density 1252 rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, pt_surface * exner_function(surface_pressure * 100.0_wp)) 1401 rho_surface = ideal_gas_law_rho(surface_pressure * 100.0_wp, & 1402 pt_surface * & 1403 exner_function(surface_pressure * 100.0_wp)) 1253 1404 ENDIF 1254 1405 … … 1273 1424 ! 1274 1425 !-- Here the seifert beheng scheme is used. Cloud concentration is assumed to 1275 !-- a constant value an qc a diagnostic value. 1426 !-- a constant value an qc a diagnostic value. 1276 1427 ELSEIF ( microphysics_seifert .AND. .NOT. microphysics_morrison ) THEN 1277 1428 CALL adjust_cloud_ij( i,j ) 1429 IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics ) & 1430 CALL ice_nucleation_ij( i,j ) 1431 IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics ) & 1432 CALL ice_deposition_ij( i,j ) 1278 1433 CALL autoconversion_ij( i,j ) 1279 1434 CALL accretion_ij( i,j ) … … 1282 1437 CALL sedimentation_rain_ij( i,j ) 1283 1438 IF ( cloud_water_sedimentation ) CALL sedimentation_cloud_ij( i,j ) 1284 1285 ! 1286 !-- Here the morrison scheme is used. No rain processes are considered and qr and nr 1439 IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics ) & 1440 CALL adjust_ice_ij ( i,j ) 1441 IF ( ice_crystal_sedimentation .AND. microphysics_ice_extension & 1442 .AND. simulated_time > start_ice_microphysics ) CALL sedimentation_ice_ij ( i,j ) 1443 ! 1444 !-- Here the morrison scheme is used. No rain processes are considered and qr and nr 1287 1445 !-- are not allocated 1288 1446 ELSEIF ( microphysics_morrison_no_rain .AND. .NOT. microphysics_seifert ) THEN 1289 1447 CALL activation_ij( i,j ) 1290 1448 CALL condensation_ij( i,j ) 1291 IF ( cloud_water_sedimentation ) CALL sedimentation_cloud_ij( i,j ) 1292 1293 ! 1294 !-- Here the full morrison scheme is used and all processes of Seifert and Beheng are 1295 !-- included 1449 IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics ) & 1450 CALL adjust_ice_ij ( i,j ) 1451 IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics ) & 1452 CALL ice_nucleation_ij( i,j ) 1453 IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics ) & 1454 CALL ice_deposition_ij( i,j ) 1455 IF ( cloud_water_sedimentation ) CALL sedimentation_cloud_ij( i,j ) 1456 1457 ! 1458 !-- Here the full morrison scheme is used and all processes of Seifert and Beheng are 1459 !-- included 1296 1460 ELSEIF ( microphysics_morrison .AND. microphysics_seifert ) THEN 1297 1461 CALL adjust_cloud_ij( i,j ) 1298 1462 CALL activation_ij( i,j ) 1299 1463 CALL condensation_ij( i,j ) 1464 IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics ) & 1465 CALL adjust_ice_ij ( i,j ) 1466 IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics ) & 1467 CALL ice_nucleation_ij( i,j ) 1468 IF ( microphysics_ice_extension .AND. simulated_time > start_ice_microphysics ) & 1469 CALL ice_deposition_ij( i,j ) 1300 1470 CALL autoconversion_ij( i,j ) 1301 1471 CALL accretion_ij( i,j ) … … 1303 1473 CALL evaporation_rain_ij( i,j ) 1304 1474 CALL sedimentation_rain_ij( i,j ) 1305 IF ( cloud_water_sedimentation ) CALL sedimentation_cloud_ij( i,j ) 1475 IF ( cloud_water_sedimentation ) CALL sedimentation_cloud_ij( i,j ) 1306 1476 1307 1477 ENDIF … … 1331 1501 IF ( microphysics_morrison ) THEN 1332 1502 CALL exchange_horiz( nc, nbgp ) 1333 CALL exchange_horiz( qc, nbgp ) 1503 CALL exchange_horiz( qc, nbgp ) 1334 1504 ENDIF 1335 1505 IF ( microphysics_seifert ) THEN … … 1337 1507 CALL exchange_horiz( nr, nbgp ) 1338 1508 ENDIF 1509 IF ( microphysics_ice_extension ) THEN 1510 CALL exchange_horiz( qi, nbgp ) 1511 CALL exchange_horiz( ni, nbgp ) 1512 ENDIF 1339 1513 CALL exchange_horiz( q, nbgp ) 1340 CALL exchange_horiz( pt, nbgp ) 1514 CALL exchange_horiz( pt, nbgp ) 1341 1515 ENDIF 1342 1516 1343 1517 1344 1518 END SUBROUTINE bcm_exchange_horiz 1345 1519 1346 1520 1347 1521 … … 1425 1599 ) & 1426 1600 * MERGE( 1.0_wp, 0.0_wp, & 1427 BTEST( wall_flags_total_0(k,j,i), 0 ) &1601 BTEST( wall_flags_total_0(k,j,i), 0 ) & 1428 1602 ) 1429 1603 IF ( qc_p(k,j,i) < 0.0_wp ) qc_p(k,j,i) = 0.0_wp … … 1517 1691 ) & 1518 1692 * MERGE( 1.0_wp, 0.0_wp, & 1519 BTEST( wall_flags_total_0(k,j,i), 0 ) &1693 BTEST( wall_flags_total_0(k,j,i), 0 ) & 1520 1694 ) 1521 1695 IF ( nc_p(k,j,i) < 0.0_wp ) nc_p(k,j,i) = 0.0_wp … … 1551 1725 1552 1726 ENDIF 1727 1728 ! 1729 !-- If required, calculate prognostic equations for ice crystal content 1730 !-- and ice crystal concentration 1731 IF ( microphysics_ice_extension ) THEN 1732 1733 CALL cpu_log( log_point(70), 'qi-equation', 'start' ) 1734 1735 ! 1736 !-- Calculate prognostic equation for ice crystal content 1737 sbt = tsc(2) 1738 IF ( scalar_advec == 'bc-scheme' ) THEN 1739 1740 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 1741 ! 1742 !-- Bott-Chlond scheme always uses Euler time step. Thus: 1743 sbt = 1.0_wp 1744 ENDIF 1745 tend = 0.0_wp 1746 CALL advec_s_bc( qi, 'qi' ) 1747 1748 ENDIF 1749 1750 ! 1751 !-- qi-tendency terms with no communication 1752 IF ( scalar_advec /= 'bc-scheme' ) THEN 1753 tend = 0.0_wp 1754 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1755 IF ( ws_scheme_sca ) THEN 1756 CALL advec_s_ws( advc_flags_s, qi, 'qi', & 1757 bc_dirichlet_l .OR. bc_radiation_l, & 1758 bc_dirichlet_n .OR. bc_radiation_n, & 1759 bc_dirichlet_r .OR. bc_radiation_r, & 1760 bc_dirichlet_s .OR. bc_radiation_s ) 1761 ELSE 1762 CALL advec_s_pw( qi ) 1763 ENDIF 1764 ELSE 1765 CALL advec_s_up( qi ) 1766 ENDIF 1767 ENDIF 1768 1769 CALL diffusion_s( qi, & 1770 surf_def_h(0)%qisws, surf_def_h(1)%qisws, & 1771 surf_def_h(2)%qisws, & 1772 surf_lsm_h%qisws, surf_usm_h%qisws, & 1773 surf_def_v(0)%qisws, surf_def_v(1)%qisws, & 1774 surf_def_v(2)%qisws, surf_def_v(3)%qisws, & 1775 surf_lsm_v(0)%qisws, surf_lsm_v(1)%qisws, & 1776 surf_lsm_v(2)%qisws, surf_lsm_v(3)%qisws, & 1777 surf_usm_v(0)%qisws, surf_usm_v(1)%qisws, & 1778 surf_usm_v(2)%qisws, surf_usm_v(3)%qisws ) 1779 1780 ! 1781 !-- Prognostic equation for ice crystal mixing ratio 1782 DO i = nxl, nxr 1783 DO j = nys, nyn 1784 !following directive is required to vectorize on Intel19 1785 !DIR$ IVDEP 1786 DO k = nzb+1, nzt 1787 qi_p(k,j,i) = qi(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + & 1788 tsc(3) * tqi_m(k,j,i) ) & 1789 - tsc(5) * rdf_sc(k) * & 1790 qi(k,j,i) & 1791 ) & 1792 * MERGE( 1.0_wp, 0.0_wp, & 1793 BTEST( wall_flags_total_0(k,j,i), 0 ) & 1794 ) 1795 IF ( qi_p(k,j,i) < 0.0_wp ) qi_p(k,j,i) = 0.0_wp 1796 ENDDO 1797 ENDDO 1798 ENDDO 1799 1800 ! 1801 !-- Calculate tendencies for the next Runge-Kutta step 1802 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1803 IF ( intermediate_timestep_count == 1 ) THEN 1804 DO i = nxl, nxr 1805 DO j = nys, nyn 1806 DO k = nzb+1, nzt 1807 tqi_m(k,j,i) = tend(k,j,i) 1808 ENDDO 1809 ENDDO 1810 ENDDO 1811 ELSEIF ( intermediate_timestep_count < & 1812 intermediate_timestep_count_max ) THEN 1813 DO i = nxl, nxr 1814 DO j = nys, nyn 1815 DO k = nzb+1, nzt 1816 tqi_m(k,j,i) = -9.5625_wp * tend(k,j,i) & 1817 + 5.3125_wp * tqi_m(k,j,i) 1818 ENDDO 1819 ENDDO 1820 ENDDO 1821 ENDIF 1822 ENDIF 1823 1824 CALL cpu_log( log_point(70), 'qi-equation', 'stop' ) 1825 1826 CALL cpu_log( log_point(69), 'ni-equation', 'start' ) 1827 ! 1828 !-- Calculate prognostic equation for ice crystal concentration 1829 sbt = tsc(2) 1830 IF ( scalar_advec == 'bc-scheme' ) THEN 1831 1832 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 1833 ! 1834 !-- Bott-Chlond scheme always uses Euler time step. Thus: 1835 sbt = 1.0_wp 1836 ENDIF 1837 tend = 0.0_wp 1838 CALL advec_s_bc( ni, 'ni' ) 1839 1840 ENDIF 1841 1842 ! 1843 !-- ni-tendency terms with no communication 1844 IF ( scalar_advec /= 'bc-scheme' ) THEN 1845 tend = 0.0_wp 1846 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1847 IF ( ws_scheme_sca ) THEN 1848 CALL advec_s_ws( advc_flags_s, ni, 'ni', & 1849 bc_dirichlet_l .OR. bc_radiation_l, & 1850 bc_dirichlet_n .OR. bc_radiation_n, & 1851 bc_dirichlet_r .OR. bc_radiation_r, & 1852 bc_dirichlet_s .OR. bc_radiation_s ) 1853 ELSE 1854 CALL advec_s_pw( ni ) 1855 ENDIF 1856 ELSE 1857 CALL advec_s_up( ni ) 1858 ENDIF 1859 ENDIF 1860 1861 CALL diffusion_s( ni, & 1862 surf_def_h(0)%nisws, surf_def_h(1)%nisws, & 1863 surf_def_h(2)%nisws, & 1864 surf_lsm_h%nisws, surf_usm_h%nisws, & 1865 surf_def_v(0)%nisws, surf_def_v(1)%nisws, & 1866 surf_def_v(2)%nisws, surf_def_v(3)%nisws, & 1867 surf_lsm_v(0)%nisws, surf_lsm_v(1)%nisws, & 1868 surf_lsm_v(2)%nisws, surf_lsm_v(3)%nisws, & 1869 surf_usm_v(0)%nisws, surf_usm_v(1)%nisws, & 1870 surf_usm_v(2)%nisws, surf_usm_v(3)%nisws ) 1871 1872 ! 1873 !-- Prognostic equation for ice crystal concentration 1874 DO i = nxl, nxr 1875 DO j = nys, nyn 1876 !following directive is required to vectorize on Intel19 1877 !DIR$ IVDEP 1878 DO k = nzb+1, nzt 1879 ni_p(k,j,i) = ni(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + & 1880 tsc(3) * tni_m(k,j,i) ) & 1881 - tsc(5) * rdf_sc(k) * & 1882 ni(k,j,i) & 1883 ) & 1884 * MERGE( 1.0_wp, 0.0_wp, & 1885 BTEST( wall_flags_total_0(k,j,i), 0 ) & 1886 ) 1887 IF ( ni_p(k,j,i) < 0.0_wp ) ni_p(k,j,i) = 0.0_wp 1888 ENDDO 1889 ENDDO 1890 ENDDO 1891 1892 ! 1893 !-- Calculate tendencies for the next Runge-Kutta step 1894 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1895 IF ( intermediate_timestep_count == 1 ) THEN 1896 DO i = nxl, nxr 1897 DO j = nys, nyn 1898 DO k = nzb+1, nzt 1899 tni_m(k,j,i) = tend(k,j,i) 1900 ENDDO 1901 ENDDO 1902 ENDDO 1903 ELSEIF ( intermediate_timestep_count < & 1904 intermediate_timestep_count_max ) THEN 1905 DO i = nxl, nxr 1906 DO j = nys, nyn 1907 DO k = nzb+1, nzt 1908 tni_m(k,j,i) = -9.5625_wp * tend(k,j,i) & 1909 + 5.3125_wp * tni_m(k,j,i) 1910 ENDDO 1911 ENDDO 1912 ENDDO 1913 ENDIF 1914 ENDIF 1915 1916 CALL cpu_log( log_point(69), 'ni-equation', 'stop' ) 1917 1918 ENDIF 1919 1553 1920 ! 1554 1921 !-- If required, calculate prognostic equations for rain water content … … 1616 1983 ) & 1617 1984 * MERGE( 1.0_wp, 0.0_wp, & 1618 BTEST( wall_flags_total_0(k,j,i), 0 ) &1985 BTEST( wall_flags_total_0(k,j,i), 0 ) & 1619 1986 ) 1620 1987 IF ( qr_p(k,j,i) < 0.0_wp ) qr_p(k,j,i) = 0.0_wp … … 1708 2075 ) & 1709 2076 * MERGE( 1.0_wp, 0.0_wp, & 1710 BTEST( wall_flags_total_0(k,j,i), 0 ) &2077 BTEST( wall_flags_total_0(k,j,i), 0 ) & 1711 2078 ) 1712 2079 IF ( nr_p(k,j,i) < 0.0_wp ) nr_p(k,j,i) = 0.0_wp … … 1751 2118 !> Control of microphysics for grid points i,j 1752 2119 !------------------------------------------------------------------------------! 1753 1754 2120 SUBROUTINE bcm_prognostic_equations_ij( i, j, i_omp_start, tn ) 1755 2121 … … 1805 2171 ) & 1806 2172 * MERGE( 1.0_wp, 0.0_wp, & 1807 BTEST( wall_flags_total_0(k,j,i), 0 )&2173 BTEST( wall_flags_total_0(k,j,i), 0 )& 1808 2174 ) 1809 2175 IF ( qc_p(k,j,i) < 0.0_wp ) qc_p(k,j,i) = 0.0_wp … … 1864 2230 ) & 1865 2231 * MERGE( 1.0_wp, 0.0_wp, & 1866 BTEST( wall_flags_total_0(k,j,i), 0 )&2232 BTEST( wall_flags_total_0(k,j,i), 0 )& 1867 2233 ) 1868 2234 IF ( nc_p(k,j,i) < 0.0_wp ) nc_p(k,j,i) = 0.0_wp … … 1885 2251 1886 2252 ENDIF 2253 2254 ! 2255 !-- If required, calculate prognostic equations for ice crystal mixing ratio 2256 !-- and ice crystal concentration 2257 IF ( microphysics_ice_extension ) THEN 2258 ! 2259 !-- Calculate prognostic equation for ice crystal mixing ratio 2260 tend(:,j,i) = 0.0_wp 2261 IF ( timestep_scheme(1:5) == 'runge' ) & 2262 THEN 2263 IF ( ws_scheme_sca ) THEN 2264 CALL advec_s_ws( advc_flags_s, i, j, qi, 'qi', flux_s_qi, & 2265 diss_s_qi, flux_l_qi, diss_l_qi, & 2266 i_omp_start, tn, & 2267 bc_dirichlet_l .OR. bc_radiation_l, & 2268 bc_dirichlet_n .OR. bc_radiation_n, & 2269 bc_dirichlet_r .OR. bc_radiation_r, & 2270 bc_dirichlet_s .OR. bc_radiation_s ) 2271 ELSE 2272 CALL advec_s_pw( i, j, qi ) 2273 ENDIF 2274 ELSE 2275 CALL advec_s_up( i, j, qi ) 2276 ENDIF 2277 CALL diffusion_s( i, j, qi, & 2278 surf_def_h(0)%qisws, surf_def_h(1)%qisws, & 2279 surf_def_h(2)%qisws, & 2280 surf_lsm_h%qisws, surf_usm_h%qisws, & 2281 surf_def_v(0)%qisws, surf_def_v(1)%qisws, & 2282 surf_def_v(2)%qisws, surf_def_v(3)%qisws, & 2283 surf_lsm_v(0)%qisws, surf_lsm_v(1)%qisws, & 2284 surf_lsm_v(2)%qisws, surf_lsm_v(3)%qisws, & 2285 surf_usm_v(0)%qisws, surf_usm_v(1)%qisws, & 2286 surf_usm_v(2)%qisws, surf_usm_v(3)%qisws ) 2287 2288 ! 2289 !-- Prognostic equation for ice crystal mixing ratio 2290 DO k = nzb+1, nzt 2291 qi_p(k,j,i) = qi(k,j,i) + ( dt_3d * & 2292 ( tsc(2) * tend(k,j,i) + & 2293 tsc(3) * tqi_m(k,j,i) )& 2294 - tsc(5) * rdf_sc(k) & 2295 * qi(k,j,i) & 2296 ) & 2297 * MERGE( 1.0_wp, 0.0_wp, & 2298 BTEST( wall_flags_total_0(k,j,i), 0 )& 2299 ) 2300 IF ( qi_p(k,j,i) < 0.0_wp ) qi_p(k,j,i) = 0.0_wp 2301 ENDDO 2302 ! 2303 !-- Calculate tendencies for the next Runge-Kutta step 2304 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2305 IF ( intermediate_timestep_count == 1 ) THEN 2306 DO k = nzb+1, nzt 2307 tqi_m(k,j,i) = tend(k,j,i) 2308 ENDDO 2309 ELSEIF ( intermediate_timestep_count < & 2310 intermediate_timestep_count_max ) THEN 2311 DO k = nzb+1, nzt 2312 tqi_m(k,j,i) = -9.5625_wp * tend(k,j,i) + & 2313 5.3125_wp * tqi_m(k,j,i) 2314 ENDDO 2315 ENDIF 2316 ENDIF 2317 2318 ! 2319 !-- Calculate prognostic equation for ice crystal concentration. 2320 tend(:,j,i) = 0.0_wp 2321 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2322 IF ( ws_scheme_sca ) THEN 2323 CALL advec_s_ws( advc_flags_s, i, j, ni, 'ni', flux_s_ni, & 2324 diss_s_ni, flux_l_ni, diss_l_ni, & 2325 i_omp_start, tn, & 2326 bc_dirichlet_l .OR. bc_radiation_l, & 2327 bc_dirichlet_n .OR. bc_radiation_n, & 2328 bc_dirichlet_r .OR. bc_radiation_r, & 2329 bc_dirichlet_s .OR. bc_radiation_s ) 2330 ELSE 2331 CALL advec_s_pw( i, j, ni ) 2332 ENDIF 2333 ELSE 2334 CALL advec_s_up( i, j, ni ) 2335 ENDIF 2336 CALL diffusion_s( i, j, ni, & 2337 surf_def_h(0)%nisws, surf_def_h(1)%nisws, & 2338 surf_def_h(2)%nisws, & 2339 surf_lsm_h%nisws, surf_usm_h%nisws, & 2340 surf_def_v(0)%nisws, surf_def_v(1)%nisws, & 2341 surf_def_v(2)%nisws, surf_def_v(3)%nisws, & 2342 surf_lsm_v(0)%nisws, surf_lsm_v(1)%nisws, & 2343 surf_lsm_v(2)%nisws, surf_lsm_v(3)%nisws, & 2344 surf_usm_v(0)%nisws, surf_usm_v(1)%nisws, & 2345 surf_usm_v(2)%nisws, surf_usm_v(3)%nisws ) 2346 2347 ! 2348 !-- Prognostic equation for ice crystal concentration 2349 DO k = nzb+1, nzt 2350 ni_p(k,j,i) = ni(k,j,i) + ( dt_3d * & 2351 ( tsc(2) * tend(k,j,i) + & 2352 tsc(3) * tni_m(k,j,i) )& 2353 - tsc(5) * rdf_sc(k) & 2354 * ni(k,j,i) & 2355 ) & 2356 * MERGE( 1.0_wp, 0.0_wp, & 2357 BTEST( wall_flags_total_0(k,j,i), 0 )& 2358 ) 2359 IF ( ni_p(k,j,i) < 0.0_wp ) ni_p(k,j,i) = 0.0_wp 2360 ENDDO 2361 ! 2362 !-- Calculate tendencies for the next Runge-Kutta step 2363 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2364 IF ( intermediate_timestep_count == 1 ) THEN 2365 DO k = nzb+1, nzt 2366 tni_m(k,j,i) = tend(k,j,i) 2367 ENDDO 2368 ELSEIF ( intermediate_timestep_count < & 2369 intermediate_timestep_count_max ) THEN 2370 DO k = nzb+1, nzt 2371 tni_m(k,j,i) = -9.5625_wp * tend(k,j,i) + & 2372 5.3125_wp * tni_m(k,j,i) 2373 ENDDO 2374 ENDIF 2375 ENDIF 2376 2377 ENDIF 2378 1887 2379 ! 1888 2380 !-- If required, calculate prognostic equations for rain water content … … 1929 2421 ) & 1930 2422 * MERGE( 1.0_wp, 0.0_wp, & 1931 BTEST( wall_flags_total_0(k,j,i), 0 )&2423 BTEST( wall_flags_total_0(k,j,i), 0 )& 1932 2424 ) 1933 2425 IF ( qr_p(k,j,i) < 0.0_wp ) qr_p(k,j,i) = 0.0_wp … … 1988 2480 ) & 1989 2481 * MERGE( 1.0_wp, 0.0_wp, & 1990 BTEST( wall_flags_total_0(k,j,i), 0 )&2482 BTEST( wall_flags_total_0(k,j,i), 0 )& 1991 2483 ) 1992 2484 IF ( nr_p(k,j,i) < 0.0_wp ) nr_p(k,j,i) = 0.0_wp … … 2038 2530 nr => nr_1; nr_p => nr_2 2039 2531 ENDIF 2532 IF ( microphysics_ice_extension ) THEN 2533 qi => qi_1; qi_p => qi_2 2534 ni => ni_1; ni_p => ni_2 2535 ENDIF 2040 2536 2041 2537 CASE ( 1 ) … … 2049 2545 nr => nr_2; nr_p => nr_1 2050 2546 ENDIF 2547 IF ( microphysics_ice_extension ) THEN 2548 qi => qi_2; qi_p => qi_1 2549 ni => ni_2; ni_p => ni_1 2550 ENDIF 2551 2051 2552 2052 2553 END SELECT … … 2093 2594 ENDIF 2094 2595 2596 IF ( microphysics_ice_extension ) THEN 2597 ! 2598 !-- Surface conditions ice crysral (Dirichlet) 2599 !-- Run loop over all non-natural and natural walls. Note, in wall-datatype 2600 !-- the k coordinate belongs to the atmospheric grid point, therefore, set 2601 !-- qr_p and nr_p at upward (k-1) and downward-facing (k+1) walls 2602 DO l = 0, 1 2603 !$OMP PARALLEL DO PRIVATE( i, j, k ) 2604 DO m = 1, bc_h(l)%ns 2605 i = bc_h(l)%i(m) 2606 j = bc_h(l)%j(m) 2607 k = bc_h(l)%k(m) 2608 qi_p(k+bc_h(l)%koff,j,i) = 0.0_wp 2609 ni_p(k+bc_h(l)%koff,j,i) = 0.0_wp 2610 ENDDO 2611 ENDDO 2612 ! 2613 !-- Top boundary condition for ice crystal (Dirichlet) 2614 qi_p(nzt+1,:,:) = 0.0_wp 2615 ni_p(nzt+1,:,:) = 0.0_wp 2616 2617 ENDIF 2618 2619 2095 2620 IF ( microphysics_seifert ) THEN 2096 2621 ! … … 2129 2654 nr_p(:,nys-1,:) = nr_p(:,nys,:) 2130 2655 ENDIF 2656 IF ( microphysics_ice_extension ) THEN 2657 qi_p(:,nys-1,:) = qi_p(:,nys,:) 2658 ni_p(:,nys-1,:) = ni_p(:,nys,:) 2659 ENDIF 2131 2660 ELSEIF ( bc_radiation_n ) THEN 2132 2661 IF ( microphysics_morrison ) THEN … … 2138 2667 nr_p(:,nyn+1,:) = nr_p(:,nyn,:) 2139 2668 ENDIF 2669 IF ( microphysics_ice_extension ) THEN 2670 qi_p(:,nyn+1,:) = qi_p(:,nyn,:) 2671 ni_p(:,nyn+1,:) = ni_p(:,nyn,:) 2672 ENDIF 2140 2673 ELSEIF ( bc_radiation_l ) THEN 2141 2674 IF ( microphysics_morrison ) THEN … … 2147 2680 nr_p(:,:,nxl-1) = nr_p(:,:,nxl) 2148 2681 ENDIF 2682 IF ( microphysics_ice_extension ) THEN 2683 qi_p(:,:,nxl-1) = qi_p(:,:,nxl) 2684 ni_p(:,:,nxl-1) = ni_p(:,:,nxl) 2685 ENDIF 2149 2686 ELSEIF ( bc_radiation_r ) THEN 2150 2687 IF ( microphysics_morrison ) THEN … … 2155 2692 qr_p(:,:,nxr+1) = qr_p(:,:,nxr) 2156 2693 nr_p(:,:,nxr+1) = nr_p(:,:,nxr) 2694 ENDIF 2695 IF ( microphysics_ice_extension ) THEN 2696 qi_p(:,:,nxr+1) = qi_p(:,:,nxr) 2697 ni_p(:,:,nxr+1) = ni_p(:,:,nxr) 2157 2698 ENDIF 2158 2699 ENDIF … … 2433 2974 ENDIF 2434 2975 to_be_resorted => nc_av 2976 ENDIF 2977 IF ( mode == 'xy' ) grid = 'zu' 2978 2979 CASE ( 'ni_xy', 'ni_xz', 'ni_yz' ) 2980 IF ( av == 0 ) THEN 2981 to_be_resorted => ni 2982 ELSE 2983 IF ( .NOT. ALLOCATED( ni_av ) ) THEN 2984 ALLOCATE( ni_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 2985 ni_av = REAL( fill_value, KIND = wp ) 2986 ENDIF 2987 to_be_resorted => ni_av 2435 2988 ENDIF 2436 2989 IF ( mode == 'xy' ) grid = 'zu' … … 2499 3052 IF ( mode == 'xy' ) grid = 'zu' 2500 3053 3054 CASE ( 'qi_xy', 'qi_xz', 'qi_yz' ) 3055 IF ( av == 0 ) THEN 3056 to_be_resorted => qi 3057 ELSE 3058 IF ( .NOT. ALLOCATED( qi_av ) ) THEN 3059 ALLOCATE( qi_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 3060 qi_av = REAL( fill_value, KIND = wp ) 3061 ENDIF 3062 to_be_resorted => qi_av 3063 ENDIF 3064 IF ( mode == 'xy' ) grid = 'zu' 3065 2501 3066 CASE ( 'ql_xy', 'ql_xz', 'ql_yz' ) 2502 3067 IF ( av == 0 ) THEN … … 2534 3099 DO k = nzb_do, nzt_do 2535 3100 local_pf(i,j,k) = MERGE( & 2536 to_be_resorted(k,j,i), 2537 REAL( fill_value, KIND = wp ), 2538 BTEST( wall_flags_total_0(k,j,i), flag_nr ) &3101 to_be_resorted(k,j,i), & 3102 REAL( fill_value, KIND = wp ), & 3103 BTEST( wall_flags_total_0(k,j,i), flag_nr ) & 2539 3104 ) 2540 3105 ENDDO … … 2595 3160 ENDIF 2596 3161 to_be_resorted => nc_av 3162 ENDIF 3163 3164 CASE ( 'ni' ) 3165 IF ( av == 0 ) THEN 3166 to_be_resorted => ni 3167 ELSE 3168 IF ( .NOT. ALLOCATED( ni_av ) ) THEN 3169 ALLOCATE( ni_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 3170 ni_av = REAL( fill_value, KIND = wp ) 3171 ENDIF 3172 to_be_resorted => ni_av 2597 3173 ENDIF 2598 3174 … … 2643 3219 ENDIF 2644 3220 3221 CASE ( 'qi' ) 3222 IF ( av == 0 ) THEN 3223 to_be_resorted => qi 3224 ELSE 3225 IF ( .NOT. ALLOCATED( qi_av ) ) THEN 3226 ALLOCATE( qi_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 3227 qi_av = REAL( fill_value, KIND = wp ) 3228 ENDIF 3229 to_be_resorted => qi_av 3230 ENDIF 3231 2645 3232 CASE ( 'ql' ) 2646 3233 IF ( av == 0 ) THEN … … 2675 3262 DO j = nys, nyn 2676 3263 DO k = nzb_do, nzt_do 2677 local_pf(i,j,k) = MERGE( 2678 to_be_resorted(k,j,i), 2679 REAL( fill_value, KIND = wp ), 2680 BTEST( wall_flags_total_0(k,j,i), flag_nr ) &3264 local_pf(i,j,k) = MERGE( & 3265 to_be_resorted(k,j,i), & 3266 REAL( fill_value, KIND = wp ), & 3267 BTEST( wall_flags_total_0(k,j,i), flag_nr ) & 2681 3268 ) 2682 3269 ENDDO … … 2750 3337 CASE ( 'curvature_solution_effects_bulk' ) 2751 3338 READ ( 13 ) curvature_solution_effects_bulk 3339 3340 CASE ( 'microphysics_ice_extension' ) 3341 READ ( 13 ) microphysics_ice_extension 3342 3343 CASE ( 'ice_crystal_sedimentation' ) 3344 READ ( 13 ) ice_crystal_sedimentation 3345 3346 CASE ( 'in_init' ) 3347 READ ( 13 ) in_init 3348 3349 CASE ( 'start_ice_microphysics' ) 3350 READ ( 13 ) start_ice_microphysics 2752 3351 2753 3352 CASE DEFAULT … … 2782 3381 CALL rrd_mpi_io( 'aerosol_bulk', aerosol_bulk ) 2783 3382 CALL rrd_mpi_io( 'curvature_solution_effects_bulk', curvature_solution_effects_bulk ) 3383 CALL rrd_mpi_io( 'start_ice_microphysics', start_ice_microphysics ) 3384 CALL rrd_mpi_io( 'microphysics_ice_extension', microphysics_ice_extension ) 3385 CALL rrd_mpi_io( 'in_init', in_init ) 3386 CALL rrd_mpi_io( 'ice_crystal_sedimentation', ice_crystal_sedimentation ) 3387 3388 2784 3389 2785 3390 END SUBROUTINE bcm_rrd_global_mpi … … 2846 3451 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 2847 3452 3453 CASE ( 'ni' ) 3454 IF ( k == 1 ) READ ( 13 ) tmp_3d 3455 ni(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 3456 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 3457 3458 CASE ( 'ni_av' ) 3459 IF ( .NOT. ALLOCATED( ni_av ) ) THEN 3460 ALLOCATE( ni_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 3461 ENDIF 3462 IF ( k == 1 ) READ ( 13 ) tmp_3d 3463 ni_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 3464 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 3465 2848 3466 CASE ( 'nr' ) 2849 3467 IF ( k == 1 ) READ ( 13 ) tmp_3d … … 2886 3504 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 2887 3505 3506 CASE ( 'qi' ) 3507 IF ( k == 1 ) READ ( 13 ) tmp_3d 3508 qi(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 3509 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 3510 3511 CASE ( 'qi_av' ) 3512 IF ( .NOT. ALLOCATED( qi_av ) ) THEN 3513 ALLOCATE( qi_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 3514 ENDIF 3515 IF ( k == 1 ) READ ( 13 ) tmp_3d 3516 qi_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = & 3517 tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp) 3518 2888 3519 CASE ( 'qc_av' ) 2889 3520 IF ( .NOT. ALLOCATED( qc_av ) ) THEN … … 2984 3615 CALL wrd_write_string( 'curvature_solution_effects_bulk' ) 2985 3616 WRITE ( 14 ) curvature_solution_effects_bulk 3617 3618 CALL wrd_write_string( 'start_ice_microphysics' ) 3619 WRITE ( 14 ) start_ice_microphysics 3620 3621 CALL wrd_write_string( 'microphysics_ice_extension' ) 3622 WRITE ( 14 ) microphysics_ice_extension 3623 3624 CALL wrd_write_string( 'in_init' ) 3625 WRITE ( 14 ) in_init 3626 3627 CALL wrd_write_string( 'ice_crystal_sedimentation' ) 3628 WRITE ( 14 ) ice_crystal_sedimentation 2986 3629 2987 3630 ELSEIF ( TRIM( restart_data_format_output ) == 'mpi' ) THEN … … 3001 3644 CALL wrd_mpi_io( 'aerosol_bulk', aerosol_bulk ) 3002 3645 CALL wrd_mpi_io( 'curvature_solution_effects_bulk', curvature_solution_effects_bulk ) 3646 CALL wrd_mpi_io( 'start_ice_microphysics', start_ice_microphysics ) 3647 CALL wrd_mpi_io( 'microphysics_ice_extension', microphysics_ice_extension ) 3648 CALL wrd_mpi_io( 'in_init', in_init ) 3649 CALL wrd_mpi_io( 'ice_crystal_sedimentation', ice_crystal_sedimentation ) 3003 3650 3004 3651 ENDIF … … 3062 3709 3063 3710 ENDIF 3711 3712 IF ( microphysics_ice_extension ) THEN 3713 3714 CALL wrd_write_string( 'ni' ) 3715 WRITE ( 14 ) ni 3716 3717 IF ( ALLOCATED( ni_av ) ) THEN 3718 CALL wrd_write_string( 'ni_av' ) 3719 WRITE ( 14 ) ni_av 3720 ENDIF 3721 3722 CALL wrd_write_string( 'qi' ) 3723 WRITE ( 14 ) qi 3724 3725 IF ( ALLOCATED( qi_av ) ) THEN 3726 CALL wrd_write_string( 'qi_av' ) 3727 WRITE ( 14 ) qi_av 3728 ENDIF 3729 3730 ENDIF 3731 3064 3732 3065 3733 IF ( microphysics_seifert ) THEN … … 3104 3772 IF ( ALLOCATED( qr_av ) ) CALL wrd_mpi_io( 'qr_av', qr_av ) 3105 3773 ENDIF 3774 IF ( microphysics_ice_extension ) THEN 3775 CALL wrd_mpi_io( 'ni', ni ) 3776 IF ( ALLOCATED( ni_av ) ) CALL wrd_mpi_io( 'ni_av', ni_av ) 3777 CALL wrd_mpi_io( 'qi', qi ) 3778 IF ( ALLOCATED( qi_av ) ) CALL wrd_mpi_io( 'qi_av', qi_av ) 3779 ENDIF 3106 3780 3107 3781 ENDIF … … 3134 3808 !-- Predetermine flag to mask topography 3135 3809 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 3136 3137 IF ( qr(k,j,i) <= eps_sb ) THEN 3138 qr(k,j,i) = 0.0_wp 3139 nr(k,j,i) = 0.0_wp 3140 ELSE 3141 IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) ) THEN 3142 nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag 3143 ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) ) THEN 3144 nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag 3810 IF ( .NOT. microphysics_morrison_no_rain ) THEN 3811 IF ( qr(k,j,i) <= eps_sb ) THEN 3812 qr(k,j,i) = 0.0_wp 3813 nr(k,j,i) = 0.0_wp 3814 ELSE 3815 IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) ) THEN 3816 nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag 3817 ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) ) THEN 3818 nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag 3819 ENDIF 3145 3820 ENDIF 3146 3821 ENDIF … … 3189 3864 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 3190 3865 3191 IF ( qr(k,j,i) <= eps_sb ) THEN 3192 qr(k,j,i) = 0.0_wp 3193 nr(k,j,i) = 0.0_wp 3194 ELSE 3195 ! 3196 !-- Adjust number of raindrops to avoid nonlinear effects in 3197 !-- sedimentation and evaporation of rain drops due to too small or 3198 !-- too big weights of rain drops (Stevens and Seifert, 2008). 3199 IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) ) THEN 3200 nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag 3201 ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) ) THEN 3202 nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag 3203 ENDIF 3204 3866 IF ( .NOT. microphysics_morrison_no_rain ) THEN 3867 IF ( qr(k,j,i) <= eps_sb ) THEN 3868 qr(k,j,i) = 0.0_wp 3869 nr(k,j,i) = 0.0_wp 3870 ELSE 3871 ! 3872 !-- Adjust number of raindrops to avoid nonlinear effects in 3873 !-- sedimentation and evaporation of rain drops due to too small or 3874 !-- too big weights of rain drops (Stevens and Seifert, 2008). 3875 IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) ) THEN 3876 nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin * flag 3877 ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) ) THEN 3878 nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax * flag 3879 ENDIF 3880 ENDIF 3205 3881 ENDIF 3206 3882 … … 3219 3895 3220 3896 END SUBROUTINE adjust_cloud_ij 3897 3898 !------------------------------------------------------------------------------! 3899 ! Description: 3900 ! ------------ 3901 !> Adjust number of ice crystal to avoid nonlinear effects in sedimentation and 3902 !> evaporation of ice crytals due to too small or too big weights 3903 !> of ice crytals (Stevens and Seifert, 2008). 3904 !------------------------------------------------------------------------------! 3905 SUBROUTINE adjust_ice 3906 3907 IMPLICIT NONE 3908 3909 INTEGER(iwp) :: i !< 3910 INTEGER(iwp) :: j !< 3911 INTEGER(iwp) :: k !< 3912 3913 REAL(wp) :: flag !< flag to indicate first grid level above 3914 3915 CALL cpu_log( log_point_s(97), 'adjust_ice', 'start' ) 3916 3917 DO i = nxl, nxr 3918 DO j = nys, nyn 3919 DO k = nzb+1, nzt 3920 ! 3921 !-- Predetermine flag to mask topography 3922 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 3923 IF ( qi(k,j,i) <= ximin ) THEN 3924 qi(k,j,i) = 0.0_wp 3925 ni(k,j,i) = 0.0_wp 3926 ELSE 3927 IF ( ni(k,j,i) * ximin > qi(k,j,i) * hyrho(k) ) THEN 3928 ni(k,j,i) = qi(k,j,i) * hyrho(k) / ximin * flag 3929 ENDIF 3930 ENDIF 3931 ENDDO 3932 ENDDO 3933 ENDDO 3934 3935 CALL cpu_log( log_point_s(97), 'adjust_ice', 'stop' ) 3936 3937 END SUBROUTINE adjust_ice 3938 3939 !------------------------------------------------------------------------------! 3940 ! Description: 3941 ! ------------ 3942 !> Adjust number of of ice crystal to avoid nonlinear effects in 3943 !> sedimentation and evaporation of ice crystals due to too small or 3944 !> too big weights of ice crytals (Stevens and Seifert, 2008). 3945 !> The same procedure is applied to cloud droplets if they are determined 3946 !> prognostically. Call for grid point i,j 3947 !------------------------------------------------------------------------------! 3948 SUBROUTINE adjust_ice_ij( i, j ) 3949 3950 IMPLICIT NONE 3951 3952 INTEGER(iwp) :: i !< 3953 INTEGER(iwp) :: j !< 3954 INTEGER(iwp) :: k !< 3955 3956 REAL(wp) :: flag !< flag to indicate first grid level above surface 3957 3958 DO k = nzb+1, nzt 3959 ! 3960 !-- Predetermine flag to mask topography 3961 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 3962 IF ( qi(k,j,i) <= ximin ) THEN 3963 qi(k,j,i) = 0.0_wp 3964 ni(k,j,i) = 0.0_wp 3965 ELSE 3966 IF ( ni(k,j,i) * ximin > qi(k,j,i) * hyrho(k) ) THEN 3967 ni(k,j,i) = qi(k,j,i) * hyrho(k) / ximin * flag 3968 ENDIF 3969 ENDIF 3970 ENDDO 3971 3972 END SUBROUTINE adjust_ice_ij 3973 3221 3974 3222 3975 !------------------------------------------------------------------------------! … … 3423 4176 END SUBROUTINE activation_ij 3424 4177 4178 !------------------------------------------------------------------------------! 4179 ! Description: 4180 ! ------------ 4181 !> Calculate ice nucleation by applying the deposition-condensation formula as 4182 !> given by Meyers et al 1992 and as described in Seifert and Beheng 2006 4183 !------------------------------------------------------------------------------! 4184 SUBROUTINE ice_nucleation 4185 4186 INTEGER(iwp) :: i !< loop index 4187 INTEGER(iwp) :: j !< loop index 4188 INTEGER(iwp) :: k !< loop index 4189 4190 LOGICAL :: isdac = .TRUE. 4191 4192 REAL(wp) :: a_m92 = -0.639_wp !< parameter for nucleation 4193 REAL(wp) :: b_m92 = 12.96_wp !< parameter for nucleation 4194 REAL(wp) :: flag !< flag to indicate first grid level 4195 REAL(wp) :: n_in !< number of ice nucleii 4196 REAL(wp) :: nucle = 0.0_wp !< nucleation rate 4197 4198 CALL cpu_log( log_point_s(93), 'ice_nucleation', 'start' ) 4199 4200 IF ( isdac ) THEN 4201 DO i = nxl, nxr 4202 DO j = nys, nyn 4203 DO k = nzb+1, nzt 4204 ! 4205 !-- Predetermine flag to mask topography 4206 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 4207 ! 4208 !-- Call calculation of supersaturation located in subroutine 4209 CALL supersaturation_ice ( i, j, k ) 4210 nucle = 0.0_wp 4211 !IF ( zu(k) >= 1500.0_wp ) CYCLE 4212 IF ( sat_ice >= 0.05_wp .OR. ql(k,j,i) >= 0.001E-3_wp ) THEN 4213 ! 4214 !-- Calculate ice nucleation 4215 nucle = MAX( ( in_init - ni(k,j,i) ) / dt_micro, 0.0_wp ) 4216 !qi(k,j,i) = qi(k,j,i) + nucle * dt_micro * ximin 4217 ni(k,j,i) = MIN( (ni(k,j,i) + nucle * dt_micro * flag), in_init) 4218 ENDIF 4219 4220 ENDDO 4221 ENDDO 4222 ENDDO 4223 ELSE 4224 DO i = nxl, nxr 4225 DO j = nys, nyn 4226 DO k = nzb+1, nzt 4227 ! 4228 !-- Predetermine flag to mask topography 4229 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 4230 ! 4231 !-- Call calculation of supersaturation located in subroutine 4232 CALL supersaturation_ice ( i, j, k ) 4233 nucle = 0.0_wp 4234 IF ( sat_ice > 0.0 ) THEN 4235 ! 4236 !-- Calculate ice nucleation 4237 n_in = in_init * EXP( a_m92 + b_m92 * sat_ice ) 4238 nucle = MAX( ( n_in - ni(k,j,i) ) / dt_micro, 0.0_wp ) 4239 ENDIF 4240 ni(k,j,i) = MIN( (ni(k,j,i) + nucle * dt_micro * flag), 1.0E10_wp ) 4241 ENDDO 4242 ENDDO 4243 ENDDO 4244 ENDIF 4245 4246 CALL cpu_log( log_point_s(93), 'ice_nucleation', 'stop' ) 4247 4248 END SUBROUTINE ice_nucleation 4249 4250 4251 !------------------------------------------------------------------------------! 4252 ! Description: 4253 ! ------------ 4254 !> Calculate ice nucleation by applying the deposition-condensation formula as 4255 !> given by Meyers et al 1992 and as described in Seifert and Beheng 2006 4256 !------------------------------------------------------------------------------! 4257 SUBROUTINE ice_nucleation_ij( i, j ) 4258 4259 INTEGER(iwp) :: i !< loop index 4260 INTEGER(iwp) :: j !< loop index 4261 INTEGER(iwp) :: k !< loop index 4262 4263 LOGICAL :: isdac = .TRUE. 4264 4265 REAL(wp) :: a_m92 = -0.639_wp !< parameter for nucleation 4266 REAL(wp) :: b_m92 = 12.96_wp !< parameter for nucleation 4267 REAL(wp) :: flag !< flag to indicate first grid level 4268 REAL(wp) :: n_in !< number of ice nucleii 4269 REAL(wp) :: nucle = 0.0_wp !< nucleation rate 4270 4271 IF ( isdac ) THEN 4272 DO k = nzb+1, nzt 4273 ! 4274 !-- Predetermine flag to mask topography 4275 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 4276 ! 4277 !-- Call calculation of supersaturation located in subroutine 4278 CALL supersaturation_ice ( i, j, k ) 4279 nucle = 0.0_wp 4280 !IF ( zu(k) >= 1500.0_wp ) CYCLE 4281 IF ( sat_ice >= 0.05_wp .OR. ql(k,j,i) >= 0.001E-3_wp ) THEN 4282 ! 4283 !-- Calculate ice nucleation 4284 nucle = MAX( ( in_init - ni(k,j,i) ) / dt_micro, 0.0_wp ) 4285 !qi(k,j,i) = qi(k,j,i) + nucle * dt_micro * ximin 4286 ni(k,j,i) = MIN( (ni(k,j,i) + nucle * dt_micro * flag), in_init) 4287 ENDIF 4288 ENDDO 4289 ELSE 4290 DO k = nzb+1, nzt 4291 ! 4292 !-- Predetermine flag to mask topography 4293 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 4294 ! 4295 !-- Call calculation of supersaturation located in subroutine 4296 CALL supersaturation_ice ( i, j, k ) 4297 nucle = 0.0_wp 4298 IF ( sat_ice > 0.0 ) THEN 4299 ! 4300 !-- Calculate ice nucleation 4301 n_in = in_init * EXP( a_m92 + b_m92 * sat_ice ) 4302 nucle = MAX( ( n_in - ni(k,j,i) ) / dt_micro, 0.0_wp ) 4303 ENDIF 4304 ni(k,j,i) = MIN( (ni(k,j,i) + nucle * dt_micro * flag), 1.0E10_wp ) 4305 ENDDO 4306 ENDIF 4307 4308 4309 END SUBROUTINE ice_nucleation_ij 3425 4310 3426 4311 !------------------------------------------------------------------------------! … … 3434 4319 IMPLICIT NONE 3435 4320 3436 INTEGER(iwp) :: i !< 3437 INTEGER(iwp) :: j !< 3438 INTEGER(iwp) :: k !< 3439 3440 REAL(wp) :: cond !< 3441 REAL(wp) :: cond_max !< 3442 REAL(wp) :: dc !< 3443 REAL(wp) :: evap !< 3444 REAL(wp) :: g_fac !< 3445 REAL(wp) :: nc_0 !< 3446 REAL(wp) :: temp !< 3447 REAL(wp) :: xc !< 3448 3449 REAL(wp) :: flag !< flag to indicate first grid level above 4321 INTEGER(iwp) :: i !< loop index 4322 INTEGER(iwp) :: j !< loop index 4323 INTEGER(iwp) :: k !< loop index 4324 4325 REAL(wp) :: cond !< condensation rate 4326 REAL(wp) :: cond_max !< maximum condensation rate 4327 REAL(wp) :: dc !< weight avageraed diameter 4328 REAL(wp) :: evap !< evaporation rate 4329 REAL(wp) :: flag !< flag to indicate first grid level above surface 4330 REAL(wp) :: g_fac !< factor 1 / Fk + Fd 4331 REAL(wp) :: nc_0 !< integral diameter 4332 REAL(wp) :: temp !< actual temperature 4333 REAL(wp) :: xc !< mean mass droplets 3450 4334 3451 4335 CALL cpu_log( log_point_s(66), 'condensation', 'start' ) … … 3461 4345 CALL supersaturation ( i, j, k ) 3462 4346 ! 3463 !-- Actual temperature: 3464 IF ( microphysics_seifert ) THEN 3465 temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) ) 3466 ELSEIF ( microphysics_morrison_no_rain ) THEN 3467 temp = t_l + lv_d_cp * qc(k,j,i) 3468 ENDIF 4347 !-- Actual temperature, t_l is calculated directly before 4348 !-- in supersaturation 4349 IF ( microphysics_ice_extension ) THEN 4350 temp = t_l + lv_d_cp * ql(k,j,i) + ls_d_cp * qi(k,j,i) 4351 ELSE 4352 temp = t_l + lv_d_cp * ql(k,j,i) 4353 ENDIF 3469 4354 3470 4355 g_fac = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * & … … 3525 4410 IMPLICIT NONE 3526 4411 3527 INTEGER(iwp) :: i !< 3528 INTEGER(iwp) :: j !< 3529 INTEGER(iwp) :: k !< 3530 3531 REAL(wp) :: cond !< 3532 REAL(wp) :: cond_max !< 3533 REAL(wp) :: dc !< 3534 REAL(wp) :: evap !< 4412 INTEGER(iwp) :: i !< loop index 4413 INTEGER(iwp) :: j !< loop index 4414 INTEGER(iwp) :: k !< loop index 4415 4416 REAL(wp) :: cond !< condensation rate 4417 REAL(wp) :: cond_max !< maximum condensation rate 4418 REAL(wp) :: dc !< weight avageraed diameter 4419 REAL(wp) :: evap !< evaporation rate 3535 4420 REAL(wp) :: flag !< flag to indicate first grid level above surface 3536 REAL(wp) :: g_fac !< 3537 REAL(wp) :: nc_0 !< 3538 REAL(wp) :: temp !< 3539 REAL(wp) :: xc !< 3540 4421 REAL(wp) :: g_fac !< factor 1 / Fk + Fd 4422 REAL(wp) :: nc_0 !< integral diameter 4423 REAL(wp) :: temp !< actual temperature 4424 REAL(wp) :: xc !< mean mass droplets 3541 4425 3542 4426 DO k = nzb+1, nzt … … 3548 4432 CALL supersaturation ( i, j, k ) 3549 4433 ! 3550 !-- Actual temperature: 3551 IF ( microphysics_seifert ) THEN 3552 temp = t_l + lv_d_cp * ( qc(k,j,i) + qr(k,j,i) ) 3553 ELSEIF ( microphysics_morrison_no_rain ) THEN 3554 temp = t_l + lv_d_cp * qc(k,j,i) 3555 ENDIF 3556 3557 g_fac = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * & 3558 l_v / ( thermal_conductivity_l * temp ) & 3559 + r_v * temp / ( diff_coeff_l * e_s ) & 4434 !-- Actual temperature, t_l is calculated directly before 4435 !-- in supersaturation 4436 IF ( microphysics_ice_extension ) THEN 4437 temp = t_l + lv_d_cp * ql(k,j,i) + ls_d_cp * qi(k,j,i) 4438 ELSE 4439 temp = t_l + lv_d_cp * ql(k,j,i) 4440 ENDIF 4441 4442 g_fac = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * & 4443 l_v / ( thermal_conductivity_l * temp ) & 4444 + r_v * temp / ( diff_coeff_l * e_s ) & 3560 4445 ) 3561 4446 ! … … 3594 4479 3595 4480 END SUBROUTINE condensation_ij 4481 4482 !------------------------------------------------------------------------------! 4483 ! Description: 4484 ! ------------ 4485 !> Calculate the growth of ice particles by water vapor deposition (after 4486 !> Seifert and Beheng, 2006). 4487 !------------------------------------------------------------------------------! 4488 SUBROUTINE ice_deposition 4489 4490 INTEGER(iwp) :: i !< loop index 4491 INTEGER(iwp) :: j !< loop index 4492 INTEGER(iwp) :: k !< loop index 4493 4494 REAL(wp) :: ac = 0.09_wp !< parameter for ice capacitance 4495 REAL(wp) :: bc = 0.33_wp !< parameter for ice capacitance 4496 REAL(wp) :: fac_gamma = 0.76_wp !< parameter to describe spectral 4497 !< distribution, here following gamma 4498 !< size distribution with µ =1/3 and nu=0 4499 REAL(wp) :: deposition_rate !< depositions rate 4500 REAL(wp) :: deposition_rate_max !< maximum deposition rate 4501 REAL(wp) :: sublimation_rate !< sublimations rate 4502 REAL(wp) :: gfac_dep !< factor 4503 REAL(wp) :: temp !< actual temperature 4504 REAL(wp) :: xi !< mean mass of ice crystal 4505 REAL(wp) :: flag !< flag to indicate first grid level above 4506 4507 CALL cpu_log( log_point_s(95), 'ice deposition', 'start' ) 4508 4509 DO i = nxl, nxr 4510 DO j = nys, nyn 4511 DO k = nzb+1, nzt 4512 ! 4513 !-- Predetermine flag to mask topography 4514 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 4515 ! 4516 !-- Call calculation of supersaturation over a plane ice surface 4517 CALL supersaturation_ice ( i, j, k ) 4518 ! 4519 !-- Actual temperature: 4520 temp = t_l + lv_d_cp * ql(k,j,i) + ls_d_cp * qi(k,j,i) 4521 4522 IF ( temp >= 273.15_wp ) CYCLE 4523 ! 4524 !-- calculating gfac_dep ( 1/ (Fk + Fd) ) see e.g. 4525 !-- Rogers and Yau, 1989 4526 gfac_dep = 1.0_wp / ( ( l_s / ( r_v * temp ) - 1.0_wp ) * & 4527 l_s / ( thermal_conductivity_l * temp ) & 4528 + r_v * temp / ( diff_coeff_l * e_si ) & 4529 ) 4530 ! 4531 !-- If there is nothing nucleated, than there is also no 4532 !-- deposition (above -38°C) 4533 IF ( ni(k,j,i) <= 0.0_wp ) CYCLE 4534 ! 4535 !-- calculate mean mass of ice crystal 4536 xi = 0.0_wp 4537 xi = MAX( ( qi(k,j,i) * hyrho(k) / ni(k,j,i)), ximin ) 4538 ! 4539 !-- Condensation needs only to be calculated in supersaturated 4540 !-- regions (regarding ice) 4541 IF ( sat_ice > 0.0_wp ) THEN 4542 ! 4543 !-- Calculate deposition rate assuming ice crystal shape as 4544 !-- prescribed in Ovchinnikov et al., 2014 and a gamma size 4545 !-- distribution according to Seifert and Beheng with to 4546 !-- µ =1/3 and nu=0 4547 deposition_rate = 4.0_wp * pi * sat_ice * gfac_dep * & 4548 fac_gamma * ac * xi**bc * ni(k,j,i) 4549 IF ( microphysics_seifert ) THEN 4550 deposition_rate_max = q(k,j,i) - & 4551 q_si - qr(k,j,i) - qi(k,j,i) 4552 ELSEIF ( microphysics_morrison_no_rain ) THEN 4553 deposition_rate_max = q(k,j,i) - & 4554 q_si - qc(k,j,i) - qi(k,j,i) 4555 ENDIF 4556 deposition_rate = MIN( deposition_rate, & 4557 deposition_rate_max / dt_micro ) 4558 4559 qi(k,j,i) = qi(k,j,i) + deposition_rate * dt_micro * flag 4560 ELSEIF ( sat_ice < 0.0_wp ) THEN 4561 sublimation_rate = 4.0_wp * pi * sat_ice * gfac_dep * & 4562 fac_gamma * ac * xi**bc * ni(k,j,i) 4563 sublimation_rate = MAX( sublimation_rate, & 4564 -qi(k,j,i) / dt_micro ) 4565 qi(k,j,i) = qi(k,j,i) + sublimation_rate * dt_micro * flag 4566 ENDIF 4567 ENDDO 4568 ENDDO 4569 ENDDO 4570 4571 CALL cpu_log( log_point_s(95), 'ice deposition', 'stop' ) 4572 4573 END SUBROUTINE ice_deposition 4574 4575 !------------------------------------------------------------------------------! 4576 ! Description: 4577 ! ------------ 4578 !> Calculate condensation rate for cloud water content (after Khairoutdinov and 4579 !> Kogan, 2000). 4580 !------------------------------------------------------------------------------! 4581 SUBROUTINE ice_deposition_ij( i, j ) 4582 4583 INTEGER(iwp) :: i !< loop index 4584 INTEGER(iwp) :: j !< loop index 4585 INTEGER(iwp) :: k !< loop index 4586 4587 REAL(wp) :: ac = 0.09_wp !< parameter for ice capacitance 4588 REAL(wp) :: bc = 0.33_wp !< parameter for ice capacitance 4589 REAL(wp) :: fac_gamma = 0.76_wp !< parameter to describe spectral 4590 !< distribution, here following gamma 4591 !< size distribution with nu=1, v=0 4592 REAL(wp) :: deposition_rate !< depositions rate 4593 REAL(wp) :: deposition_rate_max !< maximum deposition rate 4594 REAL(wp) :: sublimation_rate !< sublimations rate 4595 REAL(wp) :: gfac_dep !< factor 4596 REAL(wp) :: temp !< actual temperature 4597 REAL(wp) :: xi !< mean mass of ice crystal 4598 REAL(wp) :: flag !< flag to indicate first grid level above 4599 4600 DO k = nzb+1, nzt 4601 ! 4602 !-- Predetermine flag to mask topography 4603 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 4604 ! 4605 !-- Call calculation of supersaturation over a plane ice surface 4606 CALL supersaturation_ice ( i, j, k ) 4607 ! 4608 !-- Actual temperature: 4609 temp = t_l + lv_d_cp * ql(k,j,i) + ls_d_cp * qi(k,j,i) 4610 4611 IF ( temp >= 273.15_wp ) CYCLE 4612 ! 4613 !-- calculating gfac_dep ( 1/ (Fk + Fd) ) see e.g. 4614 !-- Rogers and Yau, 1989 4615 gfac_dep = 1.0_wp / ( ( l_s / ( r_v * temp ) - 1.0_wp ) * & 4616 l_s / ( thermal_conductivity_l * temp ) & 4617 + r_v * temp / ( diff_coeff_l * e_si ) & 4618 ) 4619 ! 4620 !-- If there is nothing nucleated, than there is also no 4621 !-- deposition (above -38°C) 4622 IF ( ni(k,j,i) <= 0.0_wp ) CYCLE 4623 ! 4624 !-- calculate mean mass of ice crystal 4625 xi = 0.0_wp 4626 xi = MAX( ( qi(k,j,i) * hyrho(k) / ni(k,j,i)), ximin ) 4627 ! 4628 !-- Condensation needs only to be calculated in supersaturated 4629 !-- regions (regarding ice) 4630 IF ( sat_ice > 0.0_wp ) THEN 4631 ! 4632 !-- Calculate deposition rate assuming ice crystal shape as 4633 !-- prescribed in Ovchinnikov et al., 2014 and a gamma size 4634 !-- distribution according to Seifert and Beheng with to 4635 !-- µ =1/3 and nu=0 4636 deposition_rate = 4.0_wp * pi * sat_ice * gfac_dep * & 4637 fac_gamma * ac * xi**bc * ni(k,j,i) 4638 IF ( microphysics_seifert ) THEN 4639 deposition_rate_max = q(k,j,i) - & 4640 q_si - qr(k,j,i) - qi(k,j,i) 4641 ELSEIF ( microphysics_morrison_no_rain ) THEN 4642 deposition_rate_max = q(k,j,i) - & 4643 q_si - qc(k,j,i) - qi(k,j,i) 4644 ENDIF 4645 deposition_rate = MIN( deposition_rate, & 4646 deposition_rate_max / dt_micro ) 4647 4648 qi(k,j,i) = qi(k,j,i) + deposition_rate * dt_micro * flag 4649 ELSEIF ( sat_ice < 0.0_wp ) THEN 4650 sublimation_rate = 4.0_wp * pi * sat_ice * gfac_dep * & 4651 fac_gamma * ac * xi**bc * ni(k,j,i) 4652 sublimation_rate = MAX( sublimation_rate, & 4653 -qi(k,j,i) / dt_micro ) 4654 qi(k,j,i) = qi(k,j,i) + sublimation_rate * dt_micro * flag 4655 ENDIF 4656 ENDDO 4657 4658 END SUBROUTINE ice_deposition_ij 3596 4659 3597 4660 … … 3823 4886 !-- Autoconversion rate (Seifert and Beheng, 2006): 3824 4887 autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / & 3825 ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 * 4888 ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 * & 3826 4889 ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) * & 3827 4890 rho_surface … … 4055 5118 ENDIF 4056 5119 4057 IF ( ( qc(k,j,i) > eps_sb ) .AND. ( qr(k,j,i) > eps_sb ) .AND. &5120 IF ( ( qc(k,j,i) > eps_sb ) .AND. ( qr(k,j,i) > eps_sb ) .AND. & 4058 5121 ( nc_accr > eps_mr ) ) THEN 4059 5122 ! … … 4082 5145 ! 4083 5146 !-- Accretion rate (Seifert and Beheng, 2006): 4084 accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac * 5147 accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac * & 4085 5148 SQRT( rho_surface * hyrho(k) ) 4086 5149 accr = MIN( accr, qc(k,j,i) / dt_micro ) … … 4089 5152 qc(k,j,i) = qc(k,j,i) - accr * dt_micro * flag 4090 5153 IF ( microphysics_morrison ) THEN 4091 nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), accr / xc * 5154 nc(k,j,i) = nc(k,j,i) - MIN( nc(k,j,i), accr / xc * & 4092 5155 hyrho(k) * dt_micro * flag & 4093 5156 ) … … 4587 5650 IF ( qc(k,j,i) > eps_sb .AND. nc(k,j,i) > eps_mr ) THEN 4588 5651 sed_nc(k) = sed_qc_const * & 4589 ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) * 5652 ( qc(k,j,i) * hyrho(k) )**( 2.0_wp / 3.0_wp ) * & 4590 5653 ( nc(k,j,i) )**( 1.0_wp / 3.0_wp ) 4591 5654 ELSE … … 4594 5657 4595 5658 sed_nc(k) = MIN( sed_nc(k), hyrho(k) * dzu(k+1) * & 4596 nc(k,j,i) / dt_micro + sed_nc(k+1) 5659 nc(k,j,i) / dt_micro + sed_nc(k+1) & 4597 5660 ) * flag 4598 5661 4599 nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) * 5662 nc(k,j,i) = nc(k,j,i) + ( sed_nc(k+1) - sed_nc(k) ) * & 4600 5663 ddzu(k+1) / hyrho(k) * dt_micro * flag 4601 5664 ENDIF … … 4608 5671 ENDIF 4609 5672 4610 sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) / 5673 sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) / & 4611 5674 dt_micro + sed_qc(k+1) & 4612 5675 ) * flag 4613 5676 4614 q(k,j,i) = q(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / 5677 q(k,j,i) = q(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / & 4615 5678 hyrho(k) * dt_micro * flag 4616 qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / 5679 qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / & 4617 5680 hyrho(k) * dt_micro * flag 4618 pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / 5681 pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / & 4619 5682 hyrho(k) * lv_d_cp * d_exner(k) * dt_micro & 4620 5683 * flag … … 4632 5695 4633 5696 END SUBROUTINE sedimentation_cloud_ij 5697 5698 !------------------------------------------------------------------------------! 5699 ! Description: 5700 ! ------------ 5701 !> Sedimentation of ice crystals 5702 !------------------------------------------------------------------------------! 5703 SUBROUTINE sedimentation_ice 5704 5705 INTEGER(iwp) :: i !< loop index 5706 INTEGER(iwp) :: j !< loop index 5707 INTEGER(iwp) :: k !< loop index 5708 5709 REAL(wp) :: flag !< flag to indicate first grid level 5710 REAL(wp) :: av = 6.39_wp !< parameter for calculating fall speed 5711 REAL(wp) :: bv = 0.1666_wp !< parameter (1/6) 5712 REAL(wp) :: xi = 0.0_wp !< mean mass of ice crystal 5713 REAL(wp) :: vi = 0.0_wp !< mean fall speed of ice crystal 5714 REAL(wp) :: factor_sed_gamma_k0 = 0.76_wp !< factor for zeroth moment and 5715 !< µ =1/3 and nu=0 5716 REAL(wp) :: factor_sed_gamma_k1 = 1.61_wp !< factor for first moment and 5717 !< µ =1/3 and nu=0 5718 5719 REAL(wp), DIMENSION(nzb:nzt+1) :: sed_ni !< sedimentation rate zeroth moment 5720 REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qi !< sedimentation rate fist moment 5721 5722 CALL cpu_log( log_point_s(96), 'sed_ice', 'start' ) 5723 5724 sed_qi(nzt+1) = 0.0_wp 5725 sed_ni(nzt+1) = 0.0_wp 5726 5727 DO i = nxl, nxr 5728 DO j = nys, nyn 5729 DO k = nzt, nzb+1, -1 5730 5731 IF ( ni(k,j,i) <= 0.0_wp ) THEN 5732 xi = 0.0_wp 5733 ELSE 5734 !-- Calculate mean mass of ice crystal 5735 xi = MAX( (hyrho(k) * qi(k,j,i) / ni(k,j,i)), ximin) 5736 ENDIF 5737 ! 5738 !-- calculate fall speed of ice crystal 5739 vi = av * xi**bv 5740 ! 5741 !-- Predetermine flag to mask topography 5742 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 5743 ! 5744 !-- Calculate sedimentation rate for each grid box, factors are 5745 !-- calculated using 5746 IF ( qi(k,j,i) > eps_sb .AND. ni(k,j,i) >= 0.0_wp ) THEN 5747 sed_qi(k) = qi(k,j,i) * vi * factor_sed_gamma_k1 * flag 5748 sed_ni(k) = ni(k,j,i) * vi * factor_sed_gamma_k0 * flag 5749 ELSE 5750 sed_qi(k) = 0.0_wp 5751 sed_ni(k) = 0.0_wp 5752 ENDIF 5753 ! 5754 !-- Calculate sedimentation: divergence of sedimentation flux 5755 sed_qi(k) = MIN( sed_qi(k), hyrho(k) * dzu(k+1) * q(k,j,i) / & 5756 dt_micro + sed_qi(k+1) & 5757 ) * flag 5758 5759 q(k,j,i) = q(k,j,i) + ( sed_qi(k+1) - sed_qi(k) ) * ddzu(k+1)& 5760 / hyrho(k) * dt_micro * flag 5761 qi(k,j,i) = qi(k,j,i) + ( sed_qi(k+1) - sed_qi(k) ) * ddzu(k+1)& 5762 / hyrho(k) * dt_micro * flag 5763 ni(k,j,i) = ni(k,j,i) + ( sed_ni(k+1) - sed_ni(k) ) * ddzu(k+1)& 5764 / hyrho(k) * dt_micro * flag 5765 5766 pt(k,j,i) = pt(k,j,i) - ( sed_qi(k+1) - sed_qi(k) ) * ddzu(k+1)& 5767 / hyrho(k) * l_s / c_p * d_exner(k) * & 5768 dt_micro * flag 5769 ! 5770 !-- Compute the precipitation rate of cloud (fog) droplets 5771 IF ( call_microphysics_at_all_substeps ) THEN 5772 prr(k,j,i) = prr(k,j,i) + sed_qi(k) / hyrho(k) * & 5773 weight_substep(intermediate_timestep_count) * & 5774 flag 5775 ELSE 5776 prr(k,j,i) = prr(k,j,i) + sed_qi(k) / hyrho(k) * flag 5777 ENDIF 5778 5779 ENDDO 5780 ENDDO 5781 ENDDO 5782 5783 CALL cpu_log( log_point_s(96), 'sed_ice', 'stop' ) 5784 5785 END SUBROUTINE sedimentation_ice 5786 5787 5788 !------------------------------------------------------------------------------! 5789 ! Description: 5790 ! ------------ 5791 !> Sedimentation of ice crystals 5792 !------------------------------------------------------------------------------! 5793 SUBROUTINE sedimentation_ice_ij( i, j ) 5794 5795 INTEGER(iwp) :: i !< 5796 INTEGER(iwp) :: j !< 5797 INTEGER(iwp) :: k !< 5798 5799 REAL(wp) :: flag !< flag to indicate first grid level 5800 REAL(wp) :: av = 6.39_wp !< parameter for calculating fall speed 5801 REAL(wp) :: bv = 0.1666_wp !< 5802 REAL(wp) :: xi = 0.0_wp !< mean mass of ice crystal 5803 REAL(wp) :: vi = 0.0_wp !< mean fall speed of ice crystal 5804 REAL(wp) :: factor_sed_gamma_k0 = 0.76_wp 5805 REAL(wp) :: factor_sed_gamma_k1 = 1.61_wp 5806 5807 REAL(wp), DIMENSION(nzb:nzt+1) :: sed_ni !< 5808 REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qi !< 5809 5810 sed_qi(nzt+1) = 0.0_wp 5811 sed_ni(nzt+1) = 0.0_wp 5812 5813 DO k = nzt, nzb+1, -1 5814 IF ( ni(k,j,i) <= 0.0_wp ) THEN 5815 xi = 0.0_wp 5816 ELSE 5817 !-- Calculate mean mass of ice crystal 5818 xi = MAX( (hyrho(k) * qi(k,j,i) / ni(k,j,i)), ximin) 5819 ENDIF 5820 ! 5821 !-- calculate fall speed of ice crystal 5822 vi = av * xi**bv 5823 ! 5824 !-- Predetermine flag to mask topography 5825 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 5826 ! 5827 !-- Calculate sedimentation rate for each grid box, factors are 5828 !-- calculated using 5829 IF ( qi(k,j,i) > eps_sb .AND. ni(k,j,i) >= 0.0_wp ) THEN 5830 sed_qi(k) = qi(k,j,i) * vi * factor_sed_gamma_k1 * flag 5831 sed_ni(k) = ni(k,j,i) * vi * factor_sed_gamma_k0 * flag 5832 ELSE 5833 sed_qi(k) = 0.0_wp 5834 sed_ni(k) = 0.0_wp 5835 ENDIF 5836 ! 5837 !-- calculate fall speed of ice crystal 5838 vi = av * xi**bv 5839 ! 5840 !-- Predetermine flag to mask topography 5841 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 5842 5843 IF ( qi(k,j,i) > eps_sb .AND. ni(k,j,i) >= 0.0_wp ) THEN 5844 sed_qi(k) = qi(k,j,i) * vi * factor_sed_gamma_k1 * flag 5845 sed_ni(k) = ni(k,j,i) * vi * factor_sed_gamma_k0 * flag 5846 ELSE 5847 sed_qi(k) = 0.0_wp 5848 sed_ni(k) = 0.0_wp 5849 ENDIF 5850 5851 sed_qi(k) = MIN( sed_qi(k), hyrho(k) * dzu(k+1) * q(k,j,i) / & 5852 dt_micro + sed_qi(k+1) & 5853 ) * flag 5854 5855 q(k,j,i) = q(k,j,i) + ( sed_qi(k+1) - sed_qi(k) ) * ddzu(k+1) / & 5856 hyrho(k) * dt_micro * flag 5857 qi(k,j,i) = qi(k,j,i) + ( sed_qi(k+1) - sed_qi(k) ) * ddzu(k+1) / & 5858 hyrho(k) * dt_micro * flag 5859 ni(k,j,i) = ni(k,j,i) + ( sed_ni(k+1) - sed_ni(k) ) * ddzu(k+1) / & 5860 hyrho(k) * dt_micro * flag 5861 pt(k,j,i) = pt(k,j,i) - ( sed_qi(k+1) - sed_qi(k) ) * ddzu(k+1) / & 5862 hyrho(k) * l_s / c_p * d_exner(k) * dt_micro & 5863 * flag 5864 ! 5865 !-- Compute the precipitation rate of cloud (fog) droplets 5866 IF ( call_microphysics_at_all_substeps ) THEN 5867 prr(k,j,i) = prr(k,j,i) + sed_qi(k) / hyrho(k) * & 5868 weight_substep(intermediate_timestep_count) * flag 5869 ELSE 5870 prr(k,j,i) = prr(k,j,i) + sed_qi(k) / hyrho(k) * flag 5871 ENDIF 5872 5873 ENDDO 5874 5875 END SUBROUTINE sedimentation_ice_ij 5876 4634 5877 4635 5878 … … 5053 6296 5054 6297 sed_nr(k) = flux / dt_micro * flag 5055 nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) / 6298 nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) / & 5056 6299 hyrho(k) * dt_micro * flag 5057 6300 ! … … 5080 6323 sed_qr(k) = flux / dt_micro * flag 5081 6324 5082 qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / 6325 qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / & 5083 6326 hyrho(k) * dt_micro * flag 5084 q(k,j,i) = q(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / 6327 q(k,j,i) = q(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / & 5085 6328 hyrho(k) * dt_micro * flag 5086 pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / 6329 pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / & 5087 6330 hyrho(k) * lv_d_cp * d_exner(k) * dt_micro & 5088 6331 * flag … … 5202 6445 !-- (see: Cuijpers + Duynkerke, 1993, JAS, 23) 5203 6446 q_s = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / ( 1.0_wp + alpha * q_s ) 6447 6448 IF ( .NOT. microphysics_ice_extension ) THEN 6449 IF ( microphysics_seifert ) THEN 6450 sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp 6451 ELSEIF ( microphysics_morrison_no_rain ) THEN 6452 sat = ( q(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp 6453 ENDIF 6454 ELSE 6455 IF ( microphysics_seifert ) THEN 6456 sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) - qi(k,j,i) ) / q_s - 1.0_wp 6457 ELSEIF ( microphysics_morrison_no_rain ) THEN 6458 sat = ( q(k,j,i) - qc(k,j,i) - qi(k,j,i) ) / q_s - 1.0_wp 6459 ENDIF 6460 ENDIF 6461 6462 END SUBROUTINE supersaturation 6463 6464 !------------------------------------------------------------------------------! 6465 ! Description: 6466 ! ------------ 6467 !> Computation of the diagnostic supersaturation sat, actual liquid water 6468 !< temperature t_l and saturation water vapor mixing ratio q_s 6469 !------------------------------------------------------------------------------! 6470 SUBROUTINE supersaturation_ice ( i, j, k ) 6471 6472 INTEGER(iwp) :: i !< running index 6473 INTEGER(iwp) :: j !< running index 6474 INTEGER(iwp) :: k !< running index 6475 6476 REAL(wp) :: e_a !< water vapor pressure 6477 REAL(wp) :: temp !< actual temperature 6478 ! 6479 !-- Actual liquid water temperature: 6480 t_l = exner(k) * pt(k,j,i) 6481 ! 6482 !-- Actual temperature: 6483 temp = t_l + lv_d_cp * ql(k,j,i) + ls_d_cp * qi(k,j,i) 6484 ! 6485 !-- Calculate water vapor saturation pressure with formular using actual 6486 !-- temperature 6487 e_si = magnus_ice( temp ) 6488 ! 6489 !-- Computation of ice saturation mixing ratio: 6490 q_si = rd_d_rv * e_si / ( hyp(k) - e_si ) 6491 ! 6492 !-- Current vapor pressure 6493 e_a = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) - qi(k,j,i) ) * hyp(k) / & 6494 ( ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) - qi(k,j,i) ) + rd_d_rv ) 5204 6495 ! 5205 6496 !-- Supersaturation: 5206 6497 !-- Not in case of microphysics_kessler or microphysics_sat_adjust 5207 6498 !-- since qr is unallocated 5208 IF ( microphysics_seifert ) THEN 5209 sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp 5210 ELSEIF ( microphysics_morrison_no_rain ) THEN 5211 sat = ( q(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp 5212 ENDIF 5213 5214 END SUBROUTINE supersaturation 6499 sat_ice = e_a / e_si - 1.0_wp 6500 6501 END SUBROUTINE supersaturation_ice 5215 6502 5216 6503 … … 5224 6511 SUBROUTINE calc_liquid_water_content 5225 6512 5226 5227 5228 IMPLICIT NONE5229 5230 6513 INTEGER(iwp) :: i !< 5231 6514 INTEGER(iwp) :: j !< … … 5234 6517 REAL(wp) :: flag !< flag to indicate first grid level above surface 5235 6518 5236 5237 6519 DO i = nxlg, nxrg 5238 6520 DO j = nysg, nyng … … 5241 6523 !-- Predetermine flag to mask topography 5242 6524 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 5243 5244 6525 ! 5245 6526 !-- Call calculation of supersaturation located 5246 6527 CALL supersaturation( i, j, k ) 5247 5248 6528 ! 5249 6529 !-- Compute the liquid water content 5250 IF ( microphysics_seifert .AND. .NOT. microphysics_morrison ) & 5251 THEN 5252 IF ( ( q(k,j,i) - q_s - qr(k,j,i) ) > 0.0_wp ) THEN 5253 qc(k,j,i) = ( q(k,j,i) - q_s - qr(k,j,i) ) * flag 5254 ql(k,j,i) = ( qc(k,j,i) + qr(k,j,i) ) * flag 6530 IF ( .NOT. microphysics_ice_extension ) THEN 6531 IF ( microphysics_seifert .AND. .NOT. & 6532 microphysics_morrison ) THEN 6533 !-- Seifert and Beheng scheme: saturation adjustment 6534 IF ( ( q(k,j,i) - q_s - qr(k,j,i) ) > 0.0_wp ) THEN 6535 qc(k,j,i) = ( q(k,j,i) - q_s - qr(k,j,i) ) * flag 6536 ql(k,j,i) = ( qc(k,j,i) + qr(k,j,i) ) * flag 6537 ELSE 6538 IF ( q(k,j,i) < qr(k,j,i) ) q(k,j,i) = qr(k,j,i) 6539 qc(k,j,i) = 0.0_wp 6540 ql(k,j,i) = qr(k,j,i) * flag 6541 ENDIF 6542 ! 6543 !-- Morrison scheme: explicit condensation (see above) 6544 ELSEIF ( microphysics_morrison .AND. & 6545 microphysics_seifert ) THEN 6546 ql(k,j,i) = qc(k,j,i) + qr(k,j,i) * flag 6547 !-- Morrison without rain: explicit condensation 6548 ELSEIF ( microphysics_morrison .AND. & 6549 .NOT. microphysics_seifert ) THEN 6550 ql(k,j,i) = qc(k,j,i) 6551 !-- Kessler and saturation adjustment scheme 5255 6552 ELSE 5256 IF ( q(k,j,i) < qr(k,j,i) ) q(k,j,i) = qr(k,j,i) 5257 qc(k,j,i) = 0.0_wp 5258 ql(k,j,i) = qr(k,j,i) * flag 6553 IF ( ( q(k,j,i) - q_s ) > 0.0_wp ) THEN 6554 qc(k,j,i) = ( q(k,j,i) - q_s ) * flag 6555 ql(k,j,i) = qc(k,j,i) * flag 6556 ELSE 6557 qc(k,j,i) = 0.0_wp 6558 ql(k,j,i) = 0.0_wp 6559 ENDIF 5259 6560 ENDIF 5260 ELSEIF ( microphysics_morrison .AND. microphysics_seifert ) THEN 5261 ql(k,j,i) = qc(k,j,i) + qr(k,j,i) * flag 5262 ELSEIF ( microphysics_morrison .AND. .NOT. microphysics_seifert ) THEN 5263 ql(k,j,i) = qc(k,j,i) 6561 !-- Calculations of liquid water content in case of mixed-phase 6562 !-- cloud microphyics 5264 6563 ELSE 5265 IF ( ( q(k,j,i) - q_s ) > 0.0_wp ) THEN 5266 qc(k,j,i) = ( q(k,j,i) - q_s ) * flag 5267 ql(k,j,i) = qc(k,j,i) * flag 6564 IF ( microphysics_seifert .AND. & 6565 .NOT. microphysics_morrison ) & 6566 THEN 6567 ! 6568 !-- Seifert and Beheng scheme: saturation adjustment 6569 IF ( ( q(k,j,i) & 6570 - q_s - qr(k,j,i) - qi(k,j,i) ) > 0.0_wp ) THEN 6571 qc(k,j,i) = ( q(k,j,i) - q_s - qr(k,j,i) - qi(k,j,i) )& 6572 * flag 6573 ql(k,j,i) = ( qc(k,j,i) + qr(k,j,i) ) * flag 6574 ELSE 6575 IF ( q(k,j,i) < ( qr(k,j,i) + qi(k,j,i) ) ) THEN 6576 q(k,j,i) = qr(k,j,i) + qi(k,j,i) 6577 ENDIF 6578 qc(k,j,i) = 0.0_wp 6579 ql(k,j,i) = qr(k,j,i) * flag 6580 ENDIF 6581 !-- Morrison scheme: explicit condensation (see above) 6582 ELSEIF ( microphysics_morrison .AND. & 6583 microphysics_seifert ) THEN 6584 ql(k,j,i) = qc(k,j,i) + qr(k,j,i) * flag 6585 !-- Morrison without rain: explicit condensation 6586 ELSEIF ( microphysics_morrison .AND. & 6587 .NOT. microphysics_seifert ) THEN 6588 ql(k,j,i) = qc(k,j,i) 6589 !-- Kessler and saturation adjustment scheme 5268 6590 ELSE 5269 qc(k,j,i) = 0.0_wp 5270 ql(k,j,i) = 0.0_wp 6591 IF ( ( q(k,j,i) - q_s ) > 0.0_wp ) THEN 6592 qc(k,j,i) = ( q(k,j,i) - q_s ) * flag 6593 ql(k,j,i) = qc(k,j,i) * flag 6594 ELSE 6595 qc(k,j,i) = 0.0_wp 6596 ql(k,j,i) = 0.0_wp 6597 ENDIF 5271 6598 ENDIF 5272 6599 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
