Changeset 1361
- Timestamp:
- Apr 16, 2014 3:17:48 PM (10 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 12 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/Makefile
r1360 r1361 20 20 # Current revisions: 21 21 # ------------------ 22 # 22 # cpulog added to microphysics 23 23 # 24 24 # Former revisions: … … 341 341 ls_forcing.o: modules.o cpulog.o mod_kinds.o 342 342 message.o: modules.o mod_kinds.o 343 microphysics.o: modules.o mod_kinds.o343 microphysics.o: modules.o cpulog.o mod_kinds.o 344 344 modules.o: modules.f90 mod_kinds.o 345 345 mod_kinds.o: mod_kinds.f90 -
palm/trunk/SOURCE/advec_s_bc.f90
r1354 r1361 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! nr and qr added 23 23 ! 24 24 ! Former revisions: … … 926 926 ENDIF 927 927 928 ELSEIF ( sk_char == 'qr' ) THEN 929 930 ! 931 !-- Rain water content boundary condition at the bottom boundary: 932 !-- Dirichlet (fixed surface rain water content). 933 DO i = nxl, nxr 934 DO j = nys, nyn 935 sk_p(nzb,j,i) = sk_p(nzb+1,j,i) 936 sk_p(nzb-1,j,i) = sk_p(nzb,j,i) 937 sk_p(nzb-2,j,i) = sk_p(nzb,j,i) 938 ENDDO 939 ENDDO 940 941 ! 942 !-- Rain water content boundary condition at the top boundary: Dirichlet 943 DO i = nxl, nxr 944 DO j = nys, nyn 945 sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i) 946 sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i) 947 ENDDO 948 ENDDO 949 950 ELSEIF ( sk_char == 'nr' ) THEN 951 952 ! 953 !-- Rain drop concentration boundary condition at the bottom boundary: 954 !-- Dirichlet (fixed surface rain drop concentration). 955 DO i = nxl, nxr 956 DO j = nys, nyn 957 sk_p(nzb,j,i) = sk_p(nzb+1,j,i) 958 sk_p(nzb-1,j,i) = sk_p(nzb,j,i) 959 sk_p(nzb-2,j,i) = sk_p(nzb,j,i) 960 ENDDO 961 ENDDO 962 963 ! 964 !-- Rain drop concentration boundary condition at the top boundary: Dirichlet 965 DO i = nxl, nxr 966 DO j = nys, nyn 967 sk_p(nzt+2,j,i) = sk_p(nzt+1,j,i) 968 sk_p(nzt+3,j,i) = sk_p(nzt+1,j,i) 969 ENDDO 970 ENDDO 971 928 972 ELSEIF ( sk_char == 'e' ) THEN 929 973 -
palm/trunk/SOURCE/advec_ws.f90
r1354 r1361 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! accelerator and vector version for qr and nr added 23 23 ! 24 24 ! Former revisions: … … 2233 2233 USE statistics, & 2234 2234 ONLY: sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l, & 2235 weight_substep2235 sums_wsqrs_ws_l, sums_wsnrs_ws_l, weight_substep 2236 2236 2237 2237 IMPLICIT NONE … … 2656 2656 * weight_substep(intermediate_timestep_count) 2657 2657 ENDDO 2658 CASE ( 'qr' ) 2659 DO k = nzb, nzt 2660 sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn) & 2661 + ( flux_t(k) + diss_t(k) ) & 2662 * weight_substep(intermediate_timestep_count) 2663 ENDDO 2664 CASE ( 'nr' ) 2665 DO k = nzb, nzt 2666 sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn) & 2667 + ( flux_t(k) + diss_t(k) ) & 2668 * weight_substep(intermediate_timestep_count) 2669 ENDDO 2658 2670 2659 2671 END SELECT … … 2690 2702 ! USE statistics, & 2691 2703 ! ONLY: sums_wspts_ws_l, sums_wsqs_ws_l, sums_wssas_ws_l, & 2692 ! weight_substep2704 ! sums_wsqrs_ws_l, sums_wsnrs_ws_l, weight_substep 2693 2705 2694 2706 IMPLICIT NONE … … 2975 2987 ! CASE ( 'q' ) 2976 2988 ! sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) & 2977 ! + ( flux_t + diss_t ) & 2989 ! + ( flux_t + diss_t ) & 2990 ! * weight_substep(intermediate_timestep_count) 2991 ! CASE ( 'qr' ) 2992 ! sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn) & 2993 ! + ( flux_t + diss_t ) & 2994 ! * weight_substep(intermediate_timestep_count) 2995 ! CASE ( 'nr' ) 2996 ! sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn) & 2997 ! + ( flux_t + diss_t ) & 2978 2998 ! * weight_substep(intermediate_timestep_count) 2979 2999 ! -
palm/trunk/SOURCE/boundary_conds.f90
r1354 r1361 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Bottom and top boundary conditions of rain water content (qr) and 23 ! rain drop concentration (nr) changed to Dirichlet 23 24 ! 24 25 ! Former revisions: … … 271 272 q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1) 272 273 273 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. & 274 precipitation ) THEN 274 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. precipitation ) THEN 275 275 ! 276 !-- Surface conditions rain water ( Neumann)276 !-- Surface conditions rain water (Dirichlet) 277 277 DO i = nxlg, nxrg 278 278 DO j = nysg, nyng 279 qr_p(nzb_s_inner(j,i),j,i) = qr_p(nzb_s_inner(j,i)+1,j,i)280 nr_p(nzb_s_inner(j,i),j,i) = nr_p(nzb_s_inner(j,i)+1,j,i)281 ENDDO 282 ENDDO 283 ! 284 !-- Top boundary condition for rain water ( Neumann)285 qr_p(nzt+1,:,:) = qr_p(nzt,:,:)286 nr_p(nzt+1,:,:) = nr_p(nzt,:,:)279 qr_p(nzb_s_inner(j,i),j,i) = 0.0_wp 280 nr_p(nzb_s_inner(j,i),j,i) = 0.0_wp 281 ENDDO 282 ENDDO 283 ! 284 !-- Top boundary condition for rain water (Dirichlet) 285 qr_p(nzt+1,:,:) = 0.0_wp 286 nr_p(nzt+1,:,:) = 0.0_wp 287 287 288 288 ENDIF -
palm/trunk/SOURCE/check_parameters.f90
r1360 r1361 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 23 ! 22 ! PA0363 removed 23 ! PA0362 changed 24 ! 24 25 ! Former revisions: 25 26 ! ----------------- … … 862 863 ENDIF 863 864 864 IF ( loop_optimization /= 'cache' .AND. cloud_physics .AND. & 865 icloud_scheme == 0 ) THEN 865 IF ( .NOT. ( loop_optimization == 'cache' .OR. & 866 loop_optimization == 'vector' ) & 867 .AND. cloud_physics .AND. icloud_scheme == 0 ) THEN 866 868 message_string = 'cloud_scheme = seifert_beheng requires ' // & 867 'loop_optimization = cache'869 'loop_optimization = "cache" or "vector"' 868 870 CALL message( 'check_parameters', 'PA0362', 1, 2, 0, 6, 0 ) 869 871 ENDIF 870 871 ! IF ( cloud_physics .AND. icloud_scheme == 0 .AND. &872 ! .NOT. precipitation .AND. .NOT. drizzle ) THEN873 ! message_string = 'cloud_scheme = seifert_beheng requires ' // &874 ! 'precipitation = .TRUE. or drizzle = .TRUE.'875 ! CALL message( 'check_parameters', 'PA0363', 1, 2, 0, 6, 0 )876 ! ENDIF877 872 878 873 ! -
palm/trunk/SOURCE/init_3d_model.f90
r1360 r1361 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! tend_* removed 23 ! Bugfix: w_subs is not allocated anymore if it is already allocated 23 24 ! 24 25 ! Former revisions: … … 366 367 ALLOCATE( qc_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 367 368 #endif 368 !369 !-- 3D-tendency arrays370 ALLOCATE( tend_pt(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &371 tend_q(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )372 369 373 370 IF ( precipitation ) THEN … … 375 372 !-- 1D-arrays 376 373 ALLOCATE ( nr_1d(nzb:nzt+1), qr_1d(nzb:nzt+1) ) 377 ! 378 ! 379 !-- 3D-tendency arrays 380 ALLOCATE( tend_nr(nzb:nzt+1,nysg:nyng,nxlg:nxrg), & 381 tend_qr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 374 382 375 ! 383 376 !-- 2D-rain water content and rain drop concentration arrays … … 475 468 ! 476 469 !-- 1D-array for large scale subsidence velocity 477 ALLOCATE ( w_subs(nzb:nzt+1) ) 478 w_subs = 0.0_wp 479 470 IF ( .NOT. ALLOCATED( w_subs ) ) THEN 471 ALLOCATE ( w_subs(nzb:nzt+1) ) 472 w_subs = 0.0_wp 473 ENDIF 480 474 481 475 ! -
palm/trunk/SOURCE/init_cloud_physics.f90
r1354 r1361 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! sed_qc_const is now calculated here (2-moment microphysics) 23 23 ! 24 24 ! Former revisions: … … 73 73 USE cloud_parameters, & 74 74 ONLY: bfactor, cp, c_sedimentation, dpirho_l, dt_precipitation, & 75 hyrho, l_d_cp, l_d_r, l_d_rv, l_v, mass_of_solute,&75 hyrho, k_st, l_d_cp, l_d_r, l_d_rv, l_v, mass_of_solute, & 76 76 molecular_weight_of_solute, molecular_weight_of_water, pirho_l, & 77 pt_d_t, rho_l, r_d, r_v, s chmidt, schmidt_p_1d3, t_d_pt,&78 vanthoff, w_precipitation77 pt_d_t, rho_l, r_d, r_v, sed_qc_const, schmidt, schmidt_p_1d3, & 78 sigma_gc, t_d_pt, vanthoff, w_precipitation 79 79 80 80 USE constants, & … … 108 108 109 109 pirho_l = pi * rho_l / 6.0_wp 110 dpirho_l = 1.0_wp / pirho_l 110 dpirho_l = 1.0_wp / pirho_l 111 ! 112 !-- constant fot the sedimentation of cloud water (2-moment cloud physics) 113 sed_qc_const = k_st * ( 3.0_wp / ( 4.0_wp * pi * rho_l ) & 114 )**( 2.0_wp / 3.0_wp ) * & 115 EXP( 5.0_wp * LOG( sigma_gc )**2 ) 111 116 ! 112 117 !-- Calculate timestep according to precipitation -
palm/trunk/SOURCE/microphysics.f90
r1354 r1361 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Bugfix in sedimentation_rain: Index corrected. 23 ! Vectorized version of adjust_cloud added. 24 ! Little reformatting of the code. 23 25 ! 24 26 ! Former revisions: … … 126 128 SUBROUTINE microphysics_control 127 129 128 USE arrays_3d 129 USE cloud_parameters 130 USE control_parameters 131 USE grid_variables 132 USE indices 130 USE arrays_3d, & 131 ONLY: hyp, nr, pt, pt_init, q, qc, qr, zu 132 133 USE cloud_parameters, & 134 ONLY: cp, hyrho, nc_const, pt_d_t, r_d, t_d_pt 135 136 USE control_parameters, & 137 ONLY: call_microphysics_at_all_substeps, drizzle, dt_3d, dt_micro, & 138 g, intermediate_timestep_count, & 139 large_scale_forcing, lsf_surf, precipitation, pt_surface, & 140 rho_surface,surface_pressure 141 142 USE indices, & 143 ONLY: nzb, nzt 144 133 145 USE kinds 134 USE statistics 146 147 USE statistics, & 148 ONLY: weight_pres 135 149 136 150 IMPLICIT NONE … … 139 153 INTEGER(iwp) :: j !: 140 154 INTEGER(iwp) :: k !: 141 142 DO i = nxl, nxr143 DO j = nys, nyn144 DO k = nzb_s_inner(j,i)+1, nzt145 146 ENDDO147 ENDDO148 ENDDO149 150 END SUBROUTINE microphysics_control151 152 SUBROUTINE adjust_cloud153 154 USE arrays_3d155 USE cloud_parameters156 USE indices157 USE kinds158 159 IMPLICIT NONE160 161 INTEGER(iwp) :: i !:162 INTEGER(iwp) :: j !:163 INTEGER(iwp) :: k !:164 165 166 DO i = nxl, nxr167 DO j = nys, nyn168 DO k = nzb_s_inner(j,i)+1, nzt169 170 ENDDO171 ENDDO172 ENDDO173 174 END SUBROUTINE adjust_cloud175 176 177 SUBROUTINE autoconversion178 179 USE arrays_3d180 USE cloud_parameters181 USE control_parameters182 USE grid_variables183 USE indices184 USE kinds185 186 IMPLICIT NONE187 188 INTEGER(iwp) :: i !:189 INTEGER(iwp) :: j !:190 INTEGER(iwp) :: k !:191 192 DO i = nxl, nxr193 DO j = nys, nyn194 DO k = nzb_s_inner(j,i)+1, nzt195 196 ENDDO197 ENDDO198 ENDDO199 200 END SUBROUTINE autoconversion201 202 203 SUBROUTINE accretion204 205 USE arrays_3d206 USE cloud_parameters207 USE control_parameters208 USE indices209 USE kinds210 211 IMPLICIT NONE212 213 INTEGER(iwp) :: i !:214 INTEGER(iwp) :: j !:215 INTEGER(iwp) :: k !:216 217 DO i = nxl, nxr218 DO j = nys, nyn219 DO k = nzb_s_inner(j,i)+1, nzt220 221 ENDDO222 ENDDO223 ENDDO224 225 END SUBROUTINE accretion226 227 228 SUBROUTINE selfcollection_breakup229 230 USE arrays_3d231 USE cloud_parameters232 USE control_parameters233 USE indices234 USE kinds235 236 IMPLICIT NONE237 238 INTEGER(iwp) :: i !:239 INTEGER(iwp) :: j !:240 INTEGER(iwp) :: k !:241 242 243 DO i = nxl, nxr244 DO j = nys, nyn245 DO k = nzb_s_inner(j,i)+1, nzt246 247 ENDDO248 ENDDO249 ENDDO250 251 END SUBROUTINE selfcollection_breakup252 253 254 SUBROUTINE evaporation_rain255 256 USE arrays_3d257 USE cloud_parameters258 USE constants259 USE control_parameters260 USE indices261 USE kinds262 263 IMPLICIT NONE264 265 INTEGER(iwp) :: i !:266 INTEGER(iwp) :: j !:267 INTEGER(iwp) :: k !:268 269 DO i = nxl, nxr270 DO j = nys, nyn271 DO k = nzb_s_inner(j,i)+1, nzt272 273 ENDDO274 ENDDO275 ENDDO276 277 END SUBROUTINE evaporation_rain278 279 280 SUBROUTINE sedimentation_cloud281 282 USE arrays_3d283 USE cloud_parameters284 USE constants285 USE control_parameters286 USE indices287 USE kinds288 289 IMPLICIT NONE290 291 INTEGER(iwp) :: i !:292 INTEGER(iwp) :: j !:293 INTEGER(iwp) :: k !:294 295 DO i = nxl, nxr296 DO j = nys, nyn297 DO k = nzb_s_inner(j,i)+1, nzt298 299 ENDDO300 ENDDO301 ENDDO302 303 END SUBROUTINE sedimentation_cloud304 305 306 SUBROUTINE sedimentation_rain307 308 USE arrays_3d309 USE cloud_parameters310 USE constants311 USE control_parameters312 USE indices313 USE kinds314 USE statistics315 316 IMPLICIT NONE317 318 INTEGER(iwp) :: i !:319 INTEGER(iwp) :: j !:320 INTEGER(iwp) :: k !:321 322 DO i = nxl, nxr323 DO j = nys, nyn324 DO k = nzb_s_inner(j,i)+1, nzt325 326 ENDDO327 ENDDO328 ENDDO329 330 END SUBROUTINE sedimentation_rain331 332 333 !------------------------------------------------------------------------------!334 ! Call for grid point i,j335 !------------------------------------------------------------------------------!336 337 SUBROUTINE microphysics_control_ij( i, j )338 339 USE arrays_3d, &340 ONLY: hyp, nc_1d, nr, nr_1d, pt, pt_init, pt_1d, q, q_1d, qc, &341 qc_1d, qr, qr_1d, tend_nr, tend_pt, tend_q, tend_qr, zu342 343 USE cloud_parameters, &344 ONLY: cp, hyrho, nc_const, pt_d_t, r_d, t_d_pt345 346 USE control_parameters, &347 ONLY: drizzle, dt_3d, dt_micro, g, intermediate_timestep_count, &348 large_scale_forcing, lsf_surf, precipitation, pt_surface, &349 rho_surface,surface_pressure350 351 USE indices, &352 ONLY: nzb, nzt353 354 USE kinds355 356 USE statistics, &357 ONLY: weight_pres358 359 IMPLICIT NONE360 361 INTEGER(iwp) :: i !:362 INTEGER(iwp) :: j !:363 INTEGER(iwp) :: k !:364 155 365 156 REAL(wp) :: t_surface !: 366 157 367 IF ( large_scale_forcing .AND.lsf_surf ) THEN158 IF ( large_scale_forcing .AND. lsf_surf ) THEN 368 159 ! 369 160 !-- Calculate: … … 374 165 DO k = nzb, nzt+1 375 166 hyp(k) = surface_pressure * 100.0_wp * & 376 ( (t_surface - g/cp * zu(k)) / t_surface )**(1.0_wp/0.286_wp) 167 ( ( t_surface - g / cp * zu(k) ) / & 168 t_surface )**(1.0_wp / 0.286_wp) 377 169 pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp 378 170 t_d_pt(k) = 1.0_wp / pt_d_t(k) … … 384 176 ENDIF 385 177 386 387 dt_micro = dt_3d * weight_pres(intermediate_timestep_count) 388 ! 389 !-- Adjust unrealistic values 390 IF ( precipitation ) CALL adjust_cloud( i,j ) 391 ! 392 !-- Use 1-d arrays 393 q_1d(:) = q(:,j,i) 394 pt_1d(:) = pt(:,j,i) 395 qc_1d(:) = qc(:,j,i) 396 nc_1d(:) = nc_const 397 IF ( precipitation ) THEN 398 qr_1d(:) = qr(:,j,i) 399 nr_1d(:) = nr(:,j,i) 178 ! 179 !-- Compute length of time step 180 IF ( call_microphysics_at_all_substeps ) THEN 181 dt_micro = dt_3d * weight_pres(intermediate_timestep_count) 182 ELSE 183 dt_micro = dt_3d 400 184 ENDIF 185 401 186 ! 402 187 !-- Compute cloud physics 403 188 IF ( precipitation ) THEN 404 CALL autoconversion( i,j ) 405 CALL accretion( i,j ) 406 CALL selfcollection_breakup( i,j ) 407 CALL evaporation_rain( i,j ) 408 CALL sedimentation_rain( i,j ) 189 CALL adjust_cloud 190 CALL autoconversion 191 CALL accretion 192 CALL selfcollection_breakup 193 CALL evaporation_rain 194 CALL sedimentation_rain 409 195 ENDIF 410 196 411 IF ( drizzle ) CALL sedimentation_cloud( i,j ) 412 ! 413 !-- Derive tendencies 414 tend_q(:,j,i) = ( q_1d(:) - q(:,j,i) ) / dt_micro 415 tend_pt(:,j,i) = ( pt_1d(:) - pt(:,j,i) ) / dt_micro 416 IF ( precipitation ) THEN 417 tend_qr(:,j,i) = ( qr_1d(:) - qr(:,j,i) ) / dt_micro 418 tend_nr(:,j,i) = ( nr_1d(:) - nr(:,j,i) ) / dt_micro 419 ENDIF 420 421 END SUBROUTINE microphysics_control_ij 422 423 SUBROUTINE adjust_cloud_ij( i, j ) 197 IF ( drizzle ) CALL sedimentation_cloud 198 199 END SUBROUTINE microphysics_control 200 201 SUBROUTINE adjust_cloud 424 202 425 203 USE arrays_3d, & … … 429 207 ONLY: eps_sb, xrmin, xrmax, hyrho, k_cc, x0 430 208 209 USE cpulog, & 210 ONLY: cpu_log, log_point_s 211 431 212 USE indices, & 432 ONLY: n zb, nzb_s_inner, nzt213 ONLY: nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt 433 214 434 215 USE kinds … … 439 220 INTEGER(iwp) :: j !: 440 221 INTEGER(iwp) :: k !: 441 ! 442 !-- Adjust number of raindrops to avoid nonlinear effects in 443 !-- sedimentation and evaporation of rain drops due to too small or 444 !-- too big weights of rain drops (Stevens and Seifert, 2008). 445 !-- The same procedure is applied to cloud droplets if they are determined 446 !-- prognostically. 447 DO k = nzb_s_inner(j,i)+1, nzt 448 449 IF ( qr(k,j,i) <= eps_sb ) THEN 450 qr(k,j,i) = 0.0_wp 451 nr(k,j,i) = 0.0_wp 452 ELSE 453 ! 454 !-- Adjust number of raindrops to avoid nonlinear effects in 455 !-- sedimentation and evaporation of rain drops due to too small or 456 !-- too big weights of rain drops (Stevens and Seifert, 2008). 457 IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) ) THEN 458 nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin 459 ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) ) THEN 460 nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax 461 ENDIF 462 463 ENDIF 464 222 223 CALL cpu_log( log_point_s(54), 'adjust_cloud', 'start' ) 224 225 DO i = nxl, nxr 226 DO j = nys, nyn 227 DO k = nzb_s_inner(j,i)+1, nzt 228 IF ( qr(k,j,i) <= eps_sb ) THEN 229 qr(k,j,i) = 0.0_wp 230 nr(k,j,i) = 0.0_wp 231 ELSE 232 ! 233 !-- Adjust number of raindrops to avoid nonlinear effects in 234 !-- sedimentation and evaporation of rain drops due to too small 235 !-- or too big weights of rain drops (Stevens and Seifert, 2008). 236 IF ( nr(k,j,i) * xrmin > qr(k,j,i) * hyrho(k) ) THEN 237 nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmin 238 ELSEIF ( nr(k,j,i) * xrmax < qr(k,j,i) * hyrho(k) ) THEN 239 nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax 240 ENDIF 241 242 ENDIF 243 ENDDO 244 ENDDO 465 245 ENDDO 466 246 467 END SUBROUTINE adjust_cloud_ij 468 469 470 SUBROUTINE autoconversion_ij( i, j ) 247 CALL cpu_log( log_point_s(54), 'adjust_cloud', 'stop' ) 248 249 END SUBROUTINE adjust_cloud 250 251 252 SUBROUTINE autoconversion 471 253 472 254 USE arrays_3d, & 473 ONLY: diss, dzu, n c_1d, nr_1d, qc_1d, qr_1d255 ONLY: diss, dzu, nr, qc, qr 474 256 475 257 USE cloud_parameters, & 476 258 ONLY: a_1, a_2, a_3, b_1, b_2, b_3, beta_cc, c_1, c_2, c_3, & 477 c_const, dpirho_l, eps_sb, hyrho, k_cc, kin_vis_air, x0 259 c_const, dpirho_l, eps_sb, hyrho, k_cc, kin_vis_air, & 260 nc_const, x0 478 261 479 262 USE control_parameters, & 480 263 ONLY: dt_micro, rho_surface, turbulence 481 264 265 USE cpulog, & 266 ONLY: cpu_log, log_point_s 267 482 268 USE grid_variables, & 483 269 ONLY: dx, dy 484 270 485 271 USE indices, & 486 ONLY: n zb, nzb_s_inner, nzt272 ONLY: nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt 487 273 488 274 USE kinds … … 496 282 REAL(wp) :: alpha_cc !: 497 283 REAL(wp) :: autocon !: 498 REAL(wp) :: epsilon!:284 REAL(wp) :: dissipation !: 499 285 REAL(wp) :: k_au !: 500 286 REAL(wp) :: l_mix !: … … 509 295 REAL(wp) :: xc !: 510 296 511 k_au = k_cc / ( 20.0_wp * x0 ) 512 513 DO k = nzb_s_inner(j,i)+1, nzt 514 515 IF ( qc_1d(k) > eps_sb ) THEN 516 ! 517 !-- Intern time scale of coagulation (Seifert and Beheng, 2006): 518 !-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr_1d(k) )) 519 tau_cloud = 1.0_wp - qc_1d(k) / ( qr_1d(k) + qc_1d(k) ) 520 ! 521 !-- Universal function for autoconversion process 522 !-- (Seifert and Beheng, 2006): 523 phi_au = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**3 524 ! 525 !-- Shape parameter of gamma distribution (Geoffroy et al., 2010): 526 !-- (Use constant nu_c = 1.0_wp instead?) 527 nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp ) 528 ! 529 !-- Mean weight of cloud droplets: 530 xc = hyrho(k) * qc_1d(k) / nc_1d(k) 531 ! 532 !-- Parameterized turbulence effects on autoconversion (Seifert, 533 !-- Nuijens and Stevens, 2010) 534 IF ( turbulence ) THEN 535 ! 536 !-- Weight averaged radius of cloud droplets: 537 rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp ) 538 539 alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c ) 540 r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c ) 541 sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c ) 542 ! 543 !-- Mixing length (neglecting distance to ground and stratification) 544 l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp ) 545 ! 546 !-- Limit dissipation rate according to Seifert, Nuijens and 547 !-- Stevens (2010) 548 epsilon = MIN( 0.06_wp, diss(k,j,i) ) 549 ! 550 !-- Compute Taylor-microscale Reynolds number: 551 re_lambda = 6.0_wp / 11.0_wp * ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) * & 552 SQRT( 15.0_wp / kin_vis_air ) * epsilon**( 1.0_wp / 6.0_wp ) 553 ! 554 !-- The factor of 1.0E4 is needed to convert the dissipation rate 555 !-- from m2 s-3 to cm2 s-3. 556 k_au = k_au * ( 1.0_wp + & 557 epsilon * 1.0E4_wp * ( re_lambda * 1.0E-3_wp )**0.25_wp * & 558 ( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) / & 559 sigma_cc )**2 ) + beta_cc ) ) 560 ENDIF 561 ! 562 !-- Autoconversion rate (Seifert and Beheng, 2006): 563 autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / & 564 ( nu_c + 1.0_wp )**2.0_wp * qc_1d(k)**2.0_wp * xc**2.0_wp * & 565 ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2.0_wp ) * & 566 rho_surface 567 autocon = MIN( autocon, qc_1d(k) / dt_micro ) 568 569 qr_1d(k) = qr_1d(k) + autocon * dt_micro 570 qc_1d(k) = qc_1d(k) - autocon * dt_micro 571 nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro 572 573 ENDIF 574 297 CALL cpu_log( log_point_s(55), 'autoconversion', 'start' ) 298 299 DO i = nxl, nxr 300 DO j = nys, nyn 301 DO k = nzb_s_inner(j,i)+1, nzt 302 303 IF ( qc(k,j,i) > eps_sb ) THEN 304 305 k_au = k_cc / ( 20.0_wp * x0 ) 306 ! 307 !-- Intern time scale of coagulation (Seifert and Beheng, 2006): 308 !-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) )) 309 tau_cloud = 1.0_wp - qc(k,j,i) / ( qr(k,j,i) + qc(k,j,i) ) 310 ! 311 !-- Universal function for autoconversion process 312 !-- (Seifert and Beheng, 2006): 313 phi_au = 600.0_wp * tau_cloud**0.68_wp * & 314 ( 1.0_wp - tau_cloud**0.68_wp )**3 315 ! 316 !-- Shape parameter of gamma distribution (Geoffroy et al., 2010): 317 !-- (Use constant nu_c = 1.0_wp instead?) 318 nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc(k,j,i) - 0.28_wp ) 319 ! 320 !-- Mean weight of cloud droplets: 321 xc = hyrho(k) * qc(k,j,i) / nc_const 322 ! 323 !-- Parameterized turbulence effects on autoconversion (Seifert, 324 !-- Nuijens and Stevens, 2010) 325 IF ( turbulence ) THEN 326 ! 327 !-- Weight averaged radius of cloud droplets: 328 rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp ) 329 330 alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c ) 331 r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c ) 332 sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c ) 333 ! 334 !-- Mixing length (neglecting distance to ground and 335 !-- stratification) 336 l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp ) 337 ! 338 !-- Limit dissipation rate according to Seifert, Nuijens and 339 !-- Stevens (2010) 340 dissipation = MIN( 0.06_wp, diss(k,j,i) ) 341 ! 342 !-- Compute Taylor-microscale Reynolds number: 343 re_lambda = 6.0_wp / 11.0_wp * & 344 ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) * & 345 SQRT( 15.0_wp / kin_vis_air ) * & 346 dissipation**( 1.0_wp / 6.0_wp ) 347 ! 348 !-- The factor of 1.0E4 is needed to convert the dissipation 349 !-- rate from m2 s-3 to cm2 s-3. 350 k_au = k_au * ( 1.0_wp + & 351 dissipation * 1.0E4_wp * & 352 ( re_lambda * 1.0E-3_wp )**0.25_wp * & 353 ( alpha_cc * EXP( -1.0_wp * ( ( rc - & 354 r_cc ) / & 355 sigma_cc )**2 & 356 ) + beta_cc & 357 ) & 358 ) 359 ENDIF 360 ! 361 !-- Autoconversion rate (Seifert and Beheng, 2006): 362 autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / & 363 ( nu_c + 1.0_wp )**2 * qc(k,j,i)**2 * xc**2 * & 364 ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) * & 365 rho_surface 366 autocon = MIN( autocon, qc(k,j,i) / dt_micro ) 367 368 qr(k,j,i) = qr(k,j,i) + autocon * dt_micro 369 qc(k,j,i) = qc(k,j,i) - autocon * dt_micro 370 nr(k,j,i) = nr(k,j,i) + autocon / x0 * hyrho(k) * dt_micro 371 372 ENDIF 373 374 ENDDO 375 ENDDO 575 376 ENDDO 576 377 577 END SUBROUTINE autoconversion_ij 578 579 580 SUBROUTINE accretion_ij( i, j ) 378 CALL cpu_log( log_point_s(55), 'autoconversion', 'stop' ) 379 380 END SUBROUTINE autoconversion 381 382 383 SUBROUTINE accretion 581 384 582 385 USE arrays_3d, & 583 ONLY: diss, qc _1d, qr_1d386 ONLY: diss, qc, qr 584 387 585 388 USE cloud_parameters, & … … 589 392 ONLY: dt_micro, rho_surface, turbulence 590 393 394 USE cpulog, & 395 ONLY: cpu_log, log_point_s 396 591 397 USE indices, & 592 ONLY: n zb, nzb_s_inner, nzt398 ONLY: nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt 593 399 594 400 USE kinds … … 606 412 REAL(wp) :: xc !: 607 413 608 DO k = nzb_s_inner(j,i)+1, nzt 609 IF ( ( qc_1d(k) > eps_sb ) .AND. ( qr_1d(k) > eps_sb ) ) THEN 610 ! 611 !-- Intern time scale of coagulation (Seifert and Beheng, 2006): 612 tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) ) 613 ! 614 !-- Universal function for accretion process 615 !-- (Seifert and Beheng, 2001): 616 phi_ac = tau_cloud / ( tau_cloud + 5.0E-5_wp ) 617 phi_ac = ( phi_ac**2.0_wp )**2.0_wp 618 ! 619 !-- Parameterized turbulence effects on autoconversion (Seifert, 620 !-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to 621 !-- convert the dissipation (diss) from m2 s-3 to cm2 s-3. 622 IF ( turbulence ) THEN 623 k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * & 624 MIN( 600.0_wp, diss(k,j,i) * 1.0E4_wp )**0.25_wp ) 625 ELSE 626 k_cr = k_cr0 627 ENDIF 628 ! 629 !-- Accretion rate (Seifert and Beheng, 2006): 630 accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac * & 631 SQRT( rho_surface * hyrho(k) ) 632 accr = MIN( accr, qc_1d(k) / dt_micro ) 633 634 qr_1d(k) = qr_1d(k) + accr * dt_micro 635 qc_1d(k) = qc_1d(k) - accr * dt_micro 636 637 ENDIF 638 414 CALL cpu_log( log_point_s(56), 'accretion', 'start' ) 415 416 DO i = nxl, nxr 417 DO j = nys, nyn 418 DO k = nzb_s_inner(j,i)+1, nzt 419 420 IF ( ( qc(k,j,i) > eps_sb ) .AND. ( qr(k,j,i) > eps_sb ) ) THEN 421 ! 422 !-- Intern time scale of coagulation (Seifert and Beheng, 2006): 423 tau_cloud = 1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr(k,j,i) ) 424 ! 425 !-- Universal function for accretion process (Seifert and 426 !-- Beheng, 2001): 427 phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4 428 ! 429 !-- Parameterized turbulence effects on autoconversion (Seifert, 430 !-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to 431 !-- convert the dissipation rate (diss) from m2 s-3 to cm2 s-3. 432 IF ( turbulence ) THEN 433 k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * & 434 MIN( 600.0_wp, & 435 diss(k,j,i) * 1.0E4_wp )**0.25_wp & 436 ) 437 ELSE 438 k_cr = k_cr0 439 ENDIF 440 ! 441 !-- Accretion rate (Seifert and Beheng, 2006): 442 accr = k_cr * qc(k,j,i) * qr(k,j,i) * phi_ac * & 443 SQRT( rho_surface * hyrho(k) ) 444 accr = MIN( accr, qc(k,j,i) / dt_micro ) 445 446 qr(k,j,i) = qr(k,j,i) + accr * dt_micro 447 qc(k,j,i) = qc(k,j,i) - accr * dt_micro 448 449 ENDIF 450 451 ENDDO 452 ENDDO 639 453 ENDDO 640 454 641 END SUBROUTINE accretion_ij 642 643 644 SUBROUTINE selfcollection_breakup_ij( i, j ) 455 CALL cpu_log( log_point_s(56), 'accretion', 'stop' ) 456 457 END SUBROUTINE accretion 458 459 460 SUBROUTINE selfcollection_breakup 645 461 646 462 USE arrays_3d, & 647 ONLY: nr _1d, qr_1d463 ONLY: nr, qr 648 464 649 465 USE cloud_parameters, & … … 653 469 ONLY: dt_micro, rho_surface 654 470 471 USE cpulog, & 472 ONLY: cpu_log, log_point_s 473 655 474 USE indices, & 656 ONLY: n zb, nzb_s_inner, nzt475 ONLY: nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt 657 476 658 477 USE kinds … … 669 488 REAL(wp) :: selfcoll !: 670 489 671 DO k = nzb_s_inner(j,i)+1, nzt 672 IF ( qr_1d(k) > eps_sb ) THEN 673 ! 674 !-- Selfcollection rate (Seifert and Beheng, 2001): 675 selfcoll = k_rr * nr_1d(k) * qr_1d(k) * & 676 SQRT( hyrho(k) * rho_surface ) 677 ! 678 !-- Weight averaged diameter of rain drops: 679 dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp ) 680 ! 681 !-- Collisional breakup rate (Seifert, 2008): 682 IF ( dr >= 0.3E-3_wp ) THEN 683 phi_br = k_br * ( dr - 1.1E-3_wp ) 684 breakup = selfcoll * ( phi_br + 1.0_wp ) 685 ELSE 686 breakup = 0.0_wp 687 ENDIF 688 689 selfcoll = MAX( breakup - selfcoll, -nr_1d(k) / dt_micro ) 690 nr_1d(k) = nr_1d(k) + selfcoll * dt_micro 691 692 ENDIF 490 CALL cpu_log( log_point_s(57), 'selfcollection', 'start' ) 491 492 DO i = nxl, nxr 493 DO j = nys, nyn 494 DO k = nzb_s_inner(j,i)+1, nzt 495 IF ( qr(k,j,i) > eps_sb ) THEN 496 ! 497 !-- Selfcollection rate (Seifert and Beheng, 2001): 498 selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * & 499 SQRT( hyrho(k) * rho_surface ) 500 ! 501 !-- Weight averaged diameter of rain drops: 502 dr = ( hyrho(k) * qr(k,j,i) / & 503 nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp ) 504 ! 505 !-- Collisional breakup rate (Seifert, 2008): 506 IF ( dr >= 0.3E-3_wp ) THEN 507 phi_br = k_br * ( dr - 1.1E-3_wp ) 508 breakup = selfcoll * ( phi_br + 1.0_wp ) 509 ELSE 510 breakup = 0.0_wp 511 ENDIF 512 513 selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / dt_micro ) 514 nr(k,j,i) = nr(k,j,i) + selfcoll * dt_micro 515 516 ENDIF 517 ENDDO 518 ENDDO 693 519 ENDDO 694 520 695 END SUBROUTINE selfcollection_breakup_ij 696 697 698 SUBROUTINE evaporation_rain_ij( i, j ) 521 CALL cpu_log( log_point_s(57), 'selfcollection', 'stop' ) 522 523 END SUBROUTINE selfcollection_breakup 524 525 526 SUBROUTINE evaporation_rain 527 699 528 ! 700 529 !-- Evaporation of precipitable water. Condensation is neglected for … … 702 531 703 532 USE arrays_3d, & 704 ONLY: hyp, nr _1d, pt_1d, q_1d, qc_1d, qr_1d533 ONLY: hyp, nr, pt, q, qc, qr 705 534 706 535 USE cloud_parameters, & … … 716 545 ONLY: dt_micro 717 546 547 USE cpulog, & 548 ONLY: cpu_log, log_point_s 549 718 550 USE indices, & 719 ONLY: n zb, nzb_s_inner, nzt551 ONLY: nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt 720 552 721 553 USE kinds … … 745 577 REAL(wp) :: xr !: 746 578 747 DO k = nzb_s_inner(j,i)+1, nzt 748 IF ( qr_1d(k) > eps_sb ) THEN 749 ! 750 !-- Actual liquid water temperature: 751 t_l = t_d_pt(k) * pt_1d(k) 752 ! 753 !-- Saturation vapor pressure at t_l: 754 e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / ( t_l - 35.86_wp ) ) 755 ! 756 !-- Computation of saturation humidity: 757 q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s ) 758 alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l ) 759 q_s = q_s * ( 1.0_wp + alpha * q_1d(k) ) / ( 1.0_wp + alpha * q_s ) 760 ! 761 !-- Supersaturation: 762 sat = MIN( 0.0_wp, ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp ) 763 ! 764 !-- Actual temperature: 765 temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) ) 766 767 g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v / & 768 ( thermal_conductivity_l * temp ) + r_v * temp / & 769 ( diff_coeff_l * e_s ) ) 770 ! 771 !-- Mean weight of rain drops 772 xr = hyrho(k) * qr_1d(k) / nr_1d(k) 773 ! 774 !-- Weight averaged diameter of rain drops: 775 dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp ) 776 ! 777 !-- Compute ventilation factor and intercept parameter 778 !-- (Seifert and Beheng, 2006; Seifert, 2008): 779 IF ( ventilation_effect ) THEN 780 ! 781 !-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005; 782 !-- Stevens and Seifert, 2008): 783 mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) ) 784 ! 785 !-- Slope parameter of gamma distribution (Seifert, 2008): 786 lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * & 787 ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr 788 789 mu_r_2 = mu_r + 2.0_wp 790 mu_r_5d2 = mu_r + 2.5_wp 791 f_vent = a_vent * gamm( mu_r_2 ) * & 792 lambda_r**( -mu_r_2 ) + & 793 b_vent * schmidt_p_1d3 * & 794 SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) * & 795 lambda_r**( -mu_r_5d2 ) * & 796 ( 1.0_wp - 0.5_wp * ( b_term / a_term ) * & 797 ( lambda_r / & 798 ( c_term + lambda_r ) )**mu_r_5d2 - & 799 0.125_wp * ( b_term / a_term )**2.0_wp * & 800 ( lambda_r / & 801 ( 2.0_wp * c_term + lambda_r ) )**mu_r_5d2 - & 802 0.0625_wp * ( b_term / a_term )**3.0_wp * & 803 ( lambda_r / & 804 ( 3.0_wp * c_term + lambda_r ) )**mu_r_5d2 - & 805 0.0390625_wp * ( b_term / a_term )**4.0_wp * & 806 ( lambda_r / & 807 ( 4.0_wp * c_term + lambda_r ) )**mu_r_5d2 ) 808 nr_0 = nr_1d(k) * lambda_r**( mu_r + 1.0_wp ) / & 809 gamm( mu_r + 1.0_wp ) 810 ELSE 811 f_vent = 1.0_wp 812 nr_0 = nr_1d(k) * dr 813 ENDIF 814 ! 815 !-- Evaporation rate of rain water content (Seifert and Beheng, 2006): 816 evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / & 817 hyrho(k) 818 819 evap = MAX( evap, -qr_1d(k) / dt_micro ) 820 evap_nr = MAX( c_evap * evap / xr * hyrho(k), & 821 -nr_1d(k) / dt_micro ) 822 823 qr_1d(k) = qr_1d(k) + evap * dt_micro 824 nr_1d(k) = nr_1d(k) + evap_nr * dt_micro 825 ENDIF 826 579 CALL cpu_log( log_point_s(58), 'evaporation', 'start' ) 580 581 DO i = nxl, nxr 582 DO j = nys, nyn 583 DO k = nzb_s_inner(j,i)+1, nzt 584 IF ( qr(k,j,i) > eps_sb ) THEN 585 ! 586 !-- Actual liquid water temperature: 587 t_l = t_d_pt(k) * pt(k,j,i) 588 ! 589 !-- Saturation vapor pressure at t_l: 590 e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / & 591 ( t_l - 35.86_wp ) & 592 ) 593 ! 594 !-- Computation of saturation humidity: 595 q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s ) 596 alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l ) 597 q_s = q_s * ( 1.0_wp + alpha * q(k,j,i) ) / & 598 ( 1.0_wp + alpha * q_s ) 599 ! 600 !-- Supersaturation: 601 sat = ( q(k,j,i) - qr(k,j,i) - qc(k,j,i) ) / q_s - 1.0_wp 602 ! 603 !-- Evaporation needs only to be calculated in subsaturated regions 604 IF ( sat < 0.0_wp ) THEN 605 ! 606 !-- Actual temperature: 607 temp = t_l + l_d_cp * ( qc(k,j,i) + qr(k,j,i) ) 608 609 g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * & 610 l_v / ( thermal_conductivity_l * temp ) & 611 + r_v * temp / ( diff_coeff_l * e_s ) & 612 ) 613 ! 614 !-- Mean weight of rain drops 615 xr = hyrho(k) * qr(k,j,i) / nr(k,j,i) 616 ! 617 !-- Weight averaged diameter of rain drops: 618 dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp ) 619 ! 620 !-- Compute ventilation factor and intercept parameter 621 !-- (Seifert and Beheng, 2006; Seifert, 2008): 622 IF ( ventilation_effect ) THEN 623 ! 624 !-- Shape parameter of gamma distribution (Milbrandt and Yau, 625 !-- 2005; Stevens and Seifert, 2008): 626 mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * & 627 ( dr - 1.4E-3_wp ) ) ) 628 ! 629 !-- Slope parameter of gamma distribution (Seifert, 2008): 630 lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * & 631 ( mu_r + 1.0_wp ) & 632 )**( 1.0_wp / 3.0_wp ) / dr 633 634 mu_r_2 = mu_r + 2.0_wp 635 mu_r_5d2 = mu_r + 2.5_wp 636 637 f_vent = a_vent * gamm( mu_r_2 ) * & 638 lambda_r**( -mu_r_2 ) + b_vent * & 639 schmidt_p_1d3 * SQRT( a_term / kin_vis_air ) *& 640 gamm( mu_r_5d2 ) * lambda_r**( -mu_r_5d2 ) * & 641 ( 1.0_wp - & 642 0.5_wp * ( b_term / a_term ) * & 643 ( lambda_r / ( c_term + lambda_r ) & 644 )**mu_r_5d2 - & 645 0.125_wp * ( b_term / a_term )**2 * & 646 ( lambda_r / ( 2.0_wp * c_term + lambda_r ) & 647 )**mu_r_5d2 - & 648 0.0625_wp * ( b_term / a_term )**3 * & 649 ( lambda_r / ( 3.0_wp * c_term + lambda_r ) & 650 )**mu_r_5d2 - & 651 0.0390625_wp * ( b_term / a_term )**4 * & 652 ( lambda_r / ( 4.0_wp * c_term + lambda_r ) & 653 )**mu_r_5d2 & 654 ) 655 656 nr_0 = nr(k,j,i) * lambda_r**( mu_r + 1.0_wp ) / & 657 gamm( mu_r + 1.0_wp ) 658 ELSE 659 f_vent = 1.0_wp 660 nr_0 = nr(k,j,i) * dr 661 ENDIF 662 ! 663 !-- Evaporation rate of rain water content (Seifert and 664 !-- Beheng, 2006): 665 evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / & 666 hyrho(k) 667 evap = MAX( evap, -qr(k,j,i) / dt_micro ) 668 evap_nr = MAX( c_evap * evap / xr * hyrho(k), & 669 -nr(k,j,i) / dt_micro ) 670 671 qr(k,j,i) = qr(k,j,i) + evap * dt_micro 672 nr(k,j,i) = nr(k,j,i) + evap_nr * dt_micro 673 674 ENDIF 675 ENDIF 676 677 ENDDO 678 ENDDO 827 679 ENDDO 828 680 829 END SUBROUTINE evaporation_rain_ij 830 831 832 SUBROUTINE sedimentation_cloud_ij( i, j ) 681 CALL cpu_log( log_point_s(58), 'evaporation', 'stop' ) 682 683 END SUBROUTINE evaporation_rain 684 685 686 SUBROUTINE sedimentation_cloud 833 687 834 688 USE arrays_3d, & 835 ONLY: ddzu, dzu, nc_1d, pt_1d, q_1d, qc_1d689 ONLY: ddzu, dzu, pt, q, qc 836 690 837 691 USE cloud_parameters, & 838 ONLY: eps_sb, hyrho, k_st, l_d_cp, prr, pt_d_t, rho_l, sigma_gc692 ONLY: eps_sb, hyrho, l_d_cp, nc_const, pt_d_t, sed_qc_const 839 693 840 694 USE constants, & … … 844 698 ONLY: dt_do2d_xy, dt_micro, intermediate_timestep_count 845 699 700 USE cpulog, & 701 ONLY: cpu_log, log_point_s 702 846 703 USE indices, & 847 ONLY: n zb, nzb_s_inner, nzt704 ONLY: nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt 848 705 849 706 USE kinds … … 855 712 INTEGER(iwp) :: k !: 856 713 857 REAL(wp) :: sed_qc_const !: 858 859 860 REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc 861 862 ! 863 !-- Sedimentation of cloud droplets (Heus et al., 2010): 864 sed_qc_const = k_st * ( 3.0_wp / ( 4.0_wp * pi * rho_l ))**( 2.0_wp / 3.0_wp ) * & 865 EXP( 5.0_wp * LOG( sigma_gc )**2 ) 866 714 REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !: 715 716 CALL cpu_log( log_point_s(59), 'sed_cloud', 'start' ) 717 718 ! 719 !-- Sedimentation of cloud droplets (Ackermann et al., 2009, MWR): 867 720 sed_qc(nzt+1) = 0.0_wp 868 721 869 DO k = nzt, nzb_s_inner(j,i)+1, -1 870 IF ( qc_1d(k) > eps_sb ) THEN 871 sed_qc(k) = sed_qc_const * nc_1d(k)**( -2.0_wp / 3.0_wp ) * & 872 ( qc_1d(k) * hyrho(k) )**( 5.0_wp / 3.0_wp ) 873 ELSE 874 sed_qc(k) = 0.0_wp 875 ENDIF 876 877 sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) / & 878 dt_micro + sed_qc(k+1) ) 879 880 q_1d(k) = q_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / & 881 hyrho(k) * dt_micro 882 qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / & 883 hyrho(k) * dt_micro 884 pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / & 885 hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro 886 722 DO i = nxl, nxr 723 DO j = nys, nyn 724 DO k = nzt, nzb_s_inner(j,i)+1, -1 725 726 IF ( qc(k,j,i) > eps_sb ) THEN 727 sed_qc(k) = sed_qc_const * nc_const**( -2.0_wp / 3.0_wp ) * & 728 ( qc(k,j,i) * hyrho(k) )**( 5.0_wp / 3.0_wp ) 729 ELSE 730 sed_qc(k) = 0.0_wp 731 ENDIF 732 733 sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q(k,j,i) / & 734 dt_micro + sed_qc(k+1) & 735 ) 736 737 q(k,j,i) = q(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * & 738 ddzu(k+1) / hyrho(k) * dt_micro 739 qc(k,j,i) = qc(k,j,i) + ( sed_qc(k+1) - sed_qc(k) ) * & 740 ddzu(k+1) / hyrho(k) * dt_micro 741 pt(k,j,i) = pt(k,j,i) - ( sed_qc(k+1) - sed_qc(k) ) * & 742 ddzu(k+1) / hyrho(k) * l_d_cp * & 743 pt_d_t(k) * dt_micro 744 745 ENDDO 746 ENDDO 887 747 ENDDO 888 748 889 END SUBROUTINE sedimentation_cloud_ij 890 891 892 SUBROUTINE sedimentation_rain_ij( i, j ) 749 CALL cpu_log( log_point_s(59), 'sed_cloud', 'stop' ) 750 751 END SUBROUTINE sedimentation_cloud 752 753 754 SUBROUTINE sedimentation_rain 893 755 894 756 USE arrays_3d, & 895 ONLY: ddzu, dzu, nr _1d, pt_1d, q_1d, qr_1d757 ONLY: ddzu, dzu, nr, pt, q, qr 896 758 897 759 USE cloud_parameters, & … … 901 763 902 764 USE control_parameters, & 903 ONLY: dt_do2d_xy, dt_micro, dt_3d, intermediate_timestep_count, & 765 ONLY: call_microphysics_at_all_substeps, dt_do2d_xy, dt_micro, & 766 dt_3d, intermediate_timestep_count, & 904 767 intermediate_timestep_count_max, & 905 768 precipitation_amount_interval, time_do2d_xy 906 769 770 USE cpulog, & 771 ONLY: cpu_log, log_point_s 772 907 773 USE indices, & 908 ONLY: n zb, nzb_s_inner, nzt774 ONLY: nxl, nxr, nys, nyn, nzb, nzb_s_inner, nzt 909 775 910 776 USE kinds … … 931 797 REAL(wp) :: z_run !: 932 798 933 REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr!:934 REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr!:935 REAL(wp), DIMENSION(nzb:nzt+1) :: d_nr!:936 REAL(wp), DIMENSION(nzb:nzt+1) :: d_qr!:937 REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope!:938 REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope!:939 REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nr!:940 REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qr!:941 REAL(wp), DIMENSION(nzb:nzt+1) :: w_nr!:942 REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr!:943 944 799 REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !: 800 REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !: 801 REAL(wp), DIMENSION(nzb:nzt+1) :: d_nr !: 802 REAL(wp), DIMENSION(nzb:nzt+1) :: d_qr !: 803 REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !: 804 REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !: 805 REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nr !: 806 REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qr !: 807 REAL(wp), DIMENSION(nzb:nzt+1) :: w_nr !: 808 REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr !: 809 810 CALL cpu_log( log_point_s(60), 'sed_rain', 'start' ) 945 811 ! 946 812 !-- Computation of sedimentation flux. Implementation according to Stevens 947 !-- and Seifert (2008). 813 !-- and Seifert (2008). Code is based on UCLA-LES. 814 IF ( intermediate_timestep_count == 1 ) prr(:,:,:) = 0.0_wp 815 ! 816 !-- Compute velocities 817 DO i = nxl, nxr 818 DO j = nys, nyn 819 DO k = nzb_s_inner(j,i)+1, nzt 820 IF ( qr(k,j,i) > eps_sb ) THEN 821 ! 822 !-- Weight averaged diameter of rain drops: 823 dr = ( hyrho(k) * qr(k,j,i) / & 824 nr(k,j,i) * dpirho_l )**( 1.0_wp / 3.0_wp ) 825 ! 826 !-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005; 827 !-- Stevens and Seifert, 2008): 828 mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * & 829 ( dr - 1.4E-3_wp ) ) ) 830 ! 831 !-- Slope parameter of gamma distribution (Seifert, 2008): 832 lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * & 833 ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr 834 835 w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, & 836 a_term - b_term * ( 1.0_wp + & 837 c_term / & 838 lambda_r )**( -1.0_wp * & 839 ( mu_r + 1.0_wp ) ) & 840 ) & 841 ) 842 843 w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, & 844 a_term - b_term * ( 1.0_wp + & 845 c_term / & 846 lambda_r )**( -1.0_wp * & 847 ( mu_r + 4.0_wp ) ) & 848 ) & 849 ) 850 ELSE 851 w_nr(k) = 0.0_wp 852 w_qr(k) = 0.0_wp 853 ENDIF 854 ENDDO 855 ! 856 !-- Adjust boundary values 857 w_nr(nzb_s_inner(j,i)) = w_nr(nzb_s_inner(j,i)+1) 858 w_qr(nzb_s_inner(j,i)) = w_qr(nzb_s_inner(j,i)+1) 859 w_nr(nzt+1) = 0.0_wp 860 w_qr(nzt+1) = 0.0_wp 861 ! 862 !-- Compute Courant number 863 DO k = nzb_s_inner(j,i)+1, nzt 864 c_nr(k) = 0.25_wp * ( w_nr(k-1) + & 865 2.0_wp * w_nr(k) + w_nr(k+1) ) * & 866 dt_micro * ddzu(k) 867 c_qr(k) = 0.25_wp * ( w_qr(k-1) + & 868 2.0_wp * w_qr(k) + w_qr(k+1) ) * & 869 dt_micro * ddzu(k) 870 ENDDO 871 ! 872 !-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977): 873 IF ( limiter_sedimentation ) THEN 874 875 DO k = nzb_s_inner(j,i)+1, nzt 876 d_mean = 0.5_wp * ( qr(k+1,j,i) + qr(k-1,j,i) ) 877 d_min = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) 878 d_max = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i) 879 880 qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, & 881 2.0_wp * d_max, & 882 ABS( d_mean ) ) 883 884 d_mean = 0.5_wp * ( nr(k+1,j,i) + nr(k-1,j,i) ) 885 d_min = nr(k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) 886 d_max = MAX( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i) 887 888 nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, & 889 2.0_wp * d_max, & 890 ABS( d_mean ) ) 891 ENDDO 892 893 ELSE 894 895 nr_slope = 0.0_wp 896 qr_slope = 0.0_wp 897 898 ENDIF 899 900 sed_nr(nzt+1) = 0.0_wp 901 sed_qr(nzt+1) = 0.0_wp 902 ! 903 !-- Compute sedimentation flux 904 DO k = nzt, nzb_s_inner(j,i)+1, -1 905 ! 906 !-- Sum up all rain drop number densities which contribute to the flux 907 !-- through k-1/2 908 flux = 0.0_wp 909 z_run = 0.0_wp ! height above z(k) 910 k_run = k 911 c_run = MIN( 1.0_wp, c_nr(k) ) 912 DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt ) 913 flux = flux + hyrho(k_run) * & 914 ( nr(k_run,j,i) + nr_slope(k_run) * & 915 ( 1.0_wp - c_run ) * 0.5_wp ) * c_run * dzu(k_run) 916 z_run = z_run + dzu(k_run) 917 k_run = k_run + 1 918 c_run = MIN( 1.0_wp, c_nr(k_run) - z_run * ddzu(k_run) ) 919 ENDDO 920 ! 921 !-- It is not allowed to sediment more rain drop number density than 922 !-- available 923 flux = MIN( flux, & 924 hyrho(k) * dzu(k+1) * nr(k,j,i) + sed_nr(k+1) * & 925 dt_micro & 926 ) 927 928 sed_nr(k) = flux / dt_micro 929 nr(k,j,i) = nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) * & 930 ddzu(k+1) / hyrho(k) * dt_micro 931 ! 932 !-- Sum up all rain water content which contributes to the flux 933 !-- through k-1/2 934 flux = 0.0_wp 935 z_run = 0.0_wp ! height above z(k) 936 k_run = k 937 c_run = MIN( 1.0_wp, c_qr(k) ) 938 939 DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt ) 940 941 flux = flux + hyrho(k_run) * ( qr(k_run,j,i) + & 942 qr_slope(k_run) * ( 1.0_wp - c_run ) * & 943 0.5_wp ) * c_run * dzu(k_run) 944 z_run = z_run + dzu(k_run) 945 k_run = k_run + 1 946 c_run = MIN( 1.0_wp, c_qr(k_run) - z_run * ddzu(k_run) ) 947 948 ENDDO 949 ! 950 !-- It is not allowed to sediment more rain water content than 951 !-- available 952 flux = MIN( flux, & 953 hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) * & 954 dt_micro & 955 ) 956 957 sed_qr(k) = flux / dt_micro 958 959 qr(k,j,i) = qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * & 960 ddzu(k+1) / hyrho(k) * dt_micro 961 q(k,j,i) = q(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * & 962 ddzu(k+1) / hyrho(k) * dt_micro 963 pt(k,j,i) = pt(k,j,i) - ( sed_qr(k+1) - sed_qr(k) ) * & 964 ddzu(k+1) / hyrho(k) * l_d_cp * & 965 pt_d_t(k) * dt_micro 966 ! 967 !-- Compute the rain rate 968 IF ( call_microphysics_at_all_substeps ) THEN 969 prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * & 970 weight_substep(intermediate_timestep_count) 971 ELSE 972 prr(k,j,i) = sed_qr(k) / hyrho(k) 973 ENDIF 974 975 ENDDO 976 ENDDO 977 ENDDO 978 979 ! 980 !-- Precipitation amount 981 IF ( intermediate_timestep_count == intermediate_timestep_count_max & 982 .AND. ( dt_do2d_xy - time_do2d_xy ) < & 983 precipitation_amount_interval ) THEN 984 DO i = nxl, nxr 985 DO j = nys, nyn 986 precipitation_amount(j,i) = precipitation_amount(j,i) + & 987 prr(nzb_s_inner(j,i)+1,j,i) * & 988 hyrho(nzb_s_inner(j,i)+1) * dt_3d 989 ENDDO 990 ENDDO 991 ENDIF 992 993 CALL cpu_log( log_point_s(60), 'sed_rain', 'stop' ) 994 995 END SUBROUTINE sedimentation_rain 996 997 998 !------------------------------------------------------------------------------! 999 ! Call for grid point i,j 1000 !------------------------------------------------------------------------------! 1001 1002 SUBROUTINE microphysics_control_ij( i, j ) 1003 1004 USE arrays_3d, & 1005 ONLY: hyp, nc_1d, nr, nr_1d, pt, pt_init, pt_1d, q, q_1d, qc, & 1006 qc_1d, qr, qr_1d, zu 1007 1008 USE cloud_parameters, & 1009 ONLY: cp, hyrho, nc_const, pt_d_t, r_d, t_d_pt 1010 1011 USE control_parameters, & 1012 ONLY: call_microphysics_at_all_substeps, drizzle, dt_3d, dt_micro, & 1013 g, intermediate_timestep_count, large_scale_forcing, & 1014 lsf_surf, precipitation, pt_surface, & 1015 rho_surface,surface_pressure 1016 1017 USE indices, & 1018 ONLY: nzb, nzt 1019 1020 USE kinds 1021 1022 USE statistics, & 1023 ONLY: weight_pres 1024 1025 IMPLICIT NONE 1026 1027 INTEGER(iwp) :: i !: 1028 INTEGER(iwp) :: j !: 1029 INTEGER(iwp) :: k !: 1030 1031 REAL(wp) :: t_surface !: 1032 1033 IF ( large_scale_forcing .AND. lsf_surf ) THEN 1034 ! 1035 !-- Calculate: 1036 !-- pt / t : ratio of potential and actual temperature (pt_d_t) 1037 !-- t / pt : ratio of actual and potential temperature (t_d_pt) 1038 !-- p_0(z) : vertical profile of the hydrostatic pressure (hyp) 1039 t_surface = pt_surface * ( surface_pressure / 1000.0_wp )**0.286_wp 1040 DO k = nzb, nzt+1 1041 hyp(k) = surface_pressure * 100.0_wp * & 1042 ( ( t_surface - g / cp * zu(k) ) / t_surface )**(1.0_wp / 0.286_wp) 1043 pt_d_t(k) = ( 100000.0_wp / hyp(k) )**0.286_wp 1044 t_d_pt(k) = 1.0_wp / pt_d_t(k) 1045 hyrho(k) = hyp(k) / ( r_d * t_d_pt(k) * pt_init(k) ) 1046 ENDDO 1047 ! 1048 !-- Compute reference density 1049 rho_surface = surface_pressure * 100.0_wp / ( r_d * t_surface ) 1050 ENDIF 1051 1052 ! 1053 !-- Compute length of time step 1054 IF ( call_microphysics_at_all_substeps ) THEN 1055 dt_micro = dt_3d * weight_pres(intermediate_timestep_count) 1056 ELSE 1057 dt_micro = dt_3d 1058 ENDIF 1059 1060 ! 1061 !-- Use 1d arrays 1062 q_1d(:) = q(:,j,i) 1063 pt_1d(:) = pt(:,j,i) 1064 qc_1d(:) = qc(:,j,i) 1065 nc_1d(:) = nc_const 1066 IF ( precipitation ) THEN 1067 qr_1d(:) = qr(:,j,i) 1068 nr_1d(:) = nr(:,j,i) 1069 ENDIF 1070 1071 ! 1072 !-- Compute cloud physics 1073 IF ( precipitation ) THEN 1074 CALL adjust_cloud( i,j ) 1075 CALL autoconversion( i,j ) 1076 CALL accretion( i,j ) 1077 CALL selfcollection_breakup( i,j ) 1078 CALL evaporation_rain( i,j ) 1079 CALL sedimentation_rain( i,j ) 1080 ENDIF 1081 1082 IF ( drizzle ) CALL sedimentation_cloud( i,j ) 1083 1084 ! 1085 !-- Store results on the 3d arrays 1086 q(:,j,i) = q_1d(:) 1087 pt(:,j,i) = pt_1d(:) 1088 IF ( precipitation ) THEN 1089 qr(:,j,i) = qr_1d(:) 1090 nr(:,j,i) = nr_1d(:) 1091 ENDIF 1092 1093 END SUBROUTINE microphysics_control_ij 1094 1095 SUBROUTINE adjust_cloud_ij( i, j ) 1096 1097 USE arrays_3d, & 1098 ONLY: qr_1d, nr_1d 1099 1100 USE cloud_parameters, & 1101 ONLY: eps_sb, xrmin, xrmax, hyrho, k_cc, x0 1102 1103 USE indices, & 1104 ONLY: nzb, nzb_s_inner, nzt 1105 1106 USE kinds 1107 1108 IMPLICIT NONE 1109 1110 INTEGER(iwp) :: i !: 1111 INTEGER(iwp) :: j !: 1112 INTEGER(iwp) :: k !: 1113 ! 1114 !-- Adjust number of raindrops to avoid nonlinear effects in 1115 !-- sedimentation and evaporation of rain drops due to too small or 1116 !-- too big weights of rain drops (Stevens and Seifert, 2008). 1117 !-- The same procedure is applied to cloud droplets if they are determined 1118 !-- prognostically. 1119 DO k = nzb_s_inner(j,i)+1, nzt 1120 1121 IF ( qr_1d(k) <= eps_sb ) THEN 1122 qr_1d(k) = 0.0_wp 1123 nr_1d(k) = 0.0_wp 1124 ELSE 1125 ! 1126 !-- Adjust number of raindrops to avoid nonlinear effects in 1127 !-- sedimentation and evaporation of rain drops due to too small or 1128 !-- too big weights of rain drops (Stevens and Seifert, 2008). 1129 IF ( nr_1d(k) * xrmin > qr_1d(k) * hyrho(k) ) THEN 1130 nr_1d(k) = qr_1d(k) * hyrho(k) / xrmin 1131 ELSEIF ( nr_1d(k) * xrmax < qr_1d(k) * hyrho(k) ) THEN 1132 nr_1d(k) = qr_1d(k) * hyrho(k) / xrmax 1133 ENDIF 1134 1135 ENDIF 1136 1137 ENDDO 1138 1139 END SUBROUTINE adjust_cloud_ij 1140 1141 1142 SUBROUTINE autoconversion_ij( i, j ) 1143 1144 USE arrays_3d, & 1145 ONLY: diss, dzu, nc_1d, nr_1d, qc_1d, qr_1d 1146 1147 USE cloud_parameters, & 1148 ONLY: a_1, a_2, a_3, b_1, b_2, b_3, beta_cc, c_1, c_2, c_3, & 1149 c_const, dpirho_l, eps_sb, hyrho, k_cc, kin_vis_air, x0 1150 1151 USE control_parameters, & 1152 ONLY: dt_micro, rho_surface, turbulence 1153 1154 USE grid_variables, & 1155 ONLY: dx, dy 1156 1157 USE indices, & 1158 ONLY: nzb, nzb_s_inner, nzt 1159 1160 USE kinds 1161 1162 IMPLICIT NONE 1163 1164 INTEGER(iwp) :: i !: 1165 INTEGER(iwp) :: j !: 1166 INTEGER(iwp) :: k !: 1167 1168 REAL(wp) :: alpha_cc !: 1169 REAL(wp) :: autocon !: 1170 REAL(wp) :: dissipation !: 1171 REAL(wp) :: k_au !: 1172 REAL(wp) :: l_mix !: 1173 REAL(wp) :: nu_c !: 1174 REAL(wp) :: phi_au !: 1175 REAL(wp) :: r_cc !: 1176 REAL(wp) :: rc !: 1177 REAL(wp) :: re_lambda !: 1178 REAL(wp) :: selfcoll !: 1179 REAL(wp) :: sigma_cc !: 1180 REAL(wp) :: tau_cloud !: 1181 REAL(wp) :: xc !: 1182 1183 DO k = nzb_s_inner(j,i)+1, nzt 1184 1185 IF ( qc_1d(k) > eps_sb ) THEN 1186 1187 k_au = k_cc / ( 20.0_wp * x0 ) 1188 ! 1189 !-- Intern time scale of coagulation (Seifert and Beheng, 2006): 1190 !-- (1.0_wp - qc(k,j,i) / ( qc(k,j,i) + qr_1d(k) )) 1191 tau_cloud = 1.0_wp - qc_1d(k) / ( qr_1d(k) + qc_1d(k) ) 1192 ! 1193 !-- Universal function for autoconversion process 1194 !-- (Seifert and Beheng, 2006): 1195 phi_au = 600.0_wp * tau_cloud**0.68_wp * ( 1.0_wp - tau_cloud**0.68_wp )**3 1196 ! 1197 !-- Shape parameter of gamma distribution (Geoffroy et al., 2010): 1198 !-- (Use constant nu_c = 1.0_wp instead?) 1199 nu_c = 1.0_wp !MAX( 0.0_wp, 1580.0_wp * hyrho(k) * qc_1d(k) - 0.28_wp ) 1200 ! 1201 !-- Mean weight of cloud droplets: 1202 xc = hyrho(k) * qc_1d(k) / nc_1d(k) 1203 ! 1204 !-- Parameterized turbulence effects on autoconversion (Seifert, 1205 !-- Nuijens and Stevens, 2010) 1206 IF ( turbulence ) THEN 1207 ! 1208 !-- Weight averaged radius of cloud droplets: 1209 rc = 0.5_wp * ( xc * dpirho_l )**( 1.0_wp / 3.0_wp ) 1210 1211 alpha_cc = ( a_1 + a_2 * nu_c ) / ( 1.0_wp + a_3 * nu_c ) 1212 r_cc = ( b_1 + b_2 * nu_c ) / ( 1.0_wp + b_3 * nu_c ) 1213 sigma_cc = ( c_1 + c_2 * nu_c ) / ( 1.0_wp + c_3 * nu_c ) 1214 ! 1215 !-- Mixing length (neglecting distance to ground and stratification) 1216 l_mix = ( dx * dy * dzu(k) )**( 1.0_wp / 3.0_wp ) 1217 ! 1218 !-- Limit dissipation rate according to Seifert, Nuijens and 1219 !-- Stevens (2010) 1220 dissipation = MIN( 0.06_wp, diss(k,j,i) ) 1221 ! 1222 !-- Compute Taylor-microscale Reynolds number: 1223 re_lambda = 6.0_wp / 11.0_wp * & 1224 ( l_mix / c_const )**( 2.0_wp / 3.0_wp ) * & 1225 SQRT( 15.0_wp / kin_vis_air ) * & 1226 dissipation**( 1.0_wp / 6.0_wp ) 1227 ! 1228 !-- The factor of 1.0E4 is needed to convert the dissipation rate 1229 !-- from m2 s-3 to cm2 s-3. 1230 k_au = k_au * ( 1.0_wp + & 1231 dissipation * 1.0E4_wp * & 1232 ( re_lambda * 1.0E-3_wp )**0.25_wp * & 1233 ( alpha_cc * EXP( -1.0_wp * ( ( rc - r_cc ) / & 1234 sigma_cc )**2 & 1235 ) + beta_cc & 1236 ) & 1237 ) 1238 ENDIF 1239 ! 1240 !-- Autoconversion rate (Seifert and Beheng, 2006): 1241 autocon = k_au * ( nu_c + 2.0_wp ) * ( nu_c + 4.0_wp ) / & 1242 ( nu_c + 1.0_wp )**2 * qc_1d(k)**2 * xc**2 * & 1243 ( 1.0_wp + phi_au / ( 1.0_wp - tau_cloud )**2 ) * & 1244 rho_surface 1245 autocon = MIN( autocon, qc_1d(k) / dt_micro ) 1246 1247 qr_1d(k) = qr_1d(k) + autocon * dt_micro 1248 qc_1d(k) = qc_1d(k) - autocon * dt_micro 1249 nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro 1250 1251 ENDIF 1252 1253 ENDDO 1254 1255 END SUBROUTINE autoconversion_ij 1256 1257 1258 SUBROUTINE accretion_ij( i, j ) 1259 1260 USE arrays_3d, & 1261 ONLY: diss, qc_1d, qr_1d 1262 1263 USE cloud_parameters, & 1264 ONLY: eps_sb, hyrho, k_cr0 1265 1266 USE control_parameters, & 1267 ONLY: dt_micro, rho_surface, turbulence 1268 1269 USE indices, & 1270 ONLY: nzb, nzb_s_inner, nzt 1271 1272 USE kinds 1273 1274 IMPLICIT NONE 1275 1276 INTEGER(iwp) :: i !: 1277 INTEGER(iwp) :: j !: 1278 INTEGER(iwp) :: k !: 1279 1280 REAL(wp) :: accr !: 1281 REAL(wp) :: k_cr !: 1282 REAL(wp) :: phi_ac !: 1283 REAL(wp) :: tau_cloud !: 1284 REAL(wp) :: xc !: 1285 1286 DO k = nzb_s_inner(j,i)+1, nzt 1287 IF ( ( qc_1d(k) > eps_sb ) .AND. ( qr_1d(k) > eps_sb ) ) THEN 1288 ! 1289 !-- Intern time scale of coagulation (Seifert and Beheng, 2006): 1290 tau_cloud = 1.0_wp - qc_1d(k) / ( qc_1d(k) + qr_1d(k) ) 1291 ! 1292 !-- Universal function for accretion process 1293 !-- (Seifert and Beheng, 2001): 1294 phi_ac = ( tau_cloud / ( tau_cloud + 5.0E-5_wp ) )**4 1295 ! 1296 !-- Parameterized turbulence effects on autoconversion (Seifert, 1297 !-- Nuijens and Stevens, 2010). The factor of 1.0E4 is needed to 1298 !-- convert the dissipation rate (diss) from m2 s-3 to cm2 s-3. 1299 IF ( turbulence ) THEN 1300 k_cr = k_cr0 * ( 1.0_wp + 0.05_wp * & 1301 MIN( 600.0_wp, & 1302 diss(k,j,i) * 1.0E4_wp )**0.25_wp & 1303 ) 1304 ELSE 1305 k_cr = k_cr0 1306 ENDIF 1307 ! 1308 !-- Accretion rate (Seifert and Beheng, 2006): 1309 accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac * SQRT( rho_surface * hyrho(k) ) 1310 accr = MIN( accr, qc_1d(k) / dt_micro ) 1311 1312 qr_1d(k) = qr_1d(k) + accr * dt_micro 1313 qc_1d(k) = qc_1d(k) - accr * dt_micro 1314 1315 ENDIF 1316 1317 ENDDO 1318 1319 END SUBROUTINE accretion_ij 1320 1321 1322 SUBROUTINE selfcollection_breakup_ij( i, j ) 1323 1324 USE arrays_3d, & 1325 ONLY: nr_1d, qr_1d 1326 1327 USE cloud_parameters, & 1328 ONLY: dpirho_l, eps_sb, hyrho, k_br, k_rr 1329 1330 USE control_parameters, & 1331 ONLY: dt_micro, rho_surface 1332 1333 USE indices, & 1334 ONLY: nzb, nzb_s_inner, nzt 1335 1336 USE kinds 1337 1338 IMPLICIT NONE 1339 1340 INTEGER(iwp) :: i !: 1341 INTEGER(iwp) :: j !: 1342 INTEGER(iwp) :: k !: 1343 1344 REAL(wp) :: breakup !: 1345 REAL(wp) :: dr !: 1346 REAL(wp) :: phi_br !: 1347 REAL(wp) :: selfcoll !: 1348 1349 DO k = nzb_s_inner(j,i)+1, nzt 1350 IF ( qr_1d(k) > eps_sb ) THEN 1351 ! 1352 !-- Selfcollection rate (Seifert and Beheng, 2001): 1353 selfcoll = k_rr * nr_1d(k) * qr_1d(k) * SQRT( hyrho(k) * rho_surface ) 1354 ! 1355 !-- Weight averaged diameter of rain drops: 1356 dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0_wp / 3.0_wp ) 1357 ! 1358 !-- Collisional breakup rate (Seifert, 2008): 1359 IF ( dr >= 0.3E-3_wp ) THEN 1360 phi_br = k_br * ( dr - 1.1E-3_wp ) 1361 breakup = selfcoll * ( phi_br + 1.0_wp ) 1362 ELSE 1363 breakup = 0.0_wp 1364 ENDIF 1365 1366 selfcoll = MAX( breakup - selfcoll, -nr_1d(k) / dt_micro ) 1367 nr_1d(k) = nr_1d(k) + selfcoll * dt_micro 1368 1369 ENDIF 1370 ENDDO 1371 1372 END SUBROUTINE selfcollection_breakup_ij 1373 1374 1375 SUBROUTINE evaporation_rain_ij( i, j ) 1376 ! 1377 !-- Evaporation of precipitable water. Condensation is neglected for 1378 !-- precipitable water. 1379 1380 USE arrays_3d, & 1381 ONLY: hyp, nr_1d, pt_1d, q_1d, qc_1d, qr_1d 1382 1383 USE cloud_parameters, & 1384 ONLY: a_term, a_vent, b_term, b_vent, c_evap, c_term, diff_coeff_l,& 1385 dpirho_l, eps_sb, hyrho, kin_vis_air, k_st, l_d_cp, l_d_r, & 1386 l_v, rho_l, r_v, schmidt_p_1d3, thermal_conductivity_l, & 1387 t_d_pt, ventilation_effect 1388 1389 USE constants, & 1390 ONLY: pi 1391 1392 USE control_parameters, & 1393 ONLY: dt_micro 1394 1395 USE indices, & 1396 ONLY: nzb, nzb_s_inner, nzt 1397 1398 USE kinds 1399 1400 IMPLICIT NONE 1401 1402 INTEGER(iwp) :: i !: 1403 INTEGER(iwp) :: j !: 1404 INTEGER(iwp) :: k !: 1405 1406 REAL(wp) :: alpha !: 1407 REAL(wp) :: dr !: 1408 REAL(wp) :: e_s !: 1409 REAL(wp) :: evap !: 1410 REAL(wp) :: evap_nr !: 1411 REAL(wp) :: f_vent !: 1412 REAL(wp) :: g_evap !: 1413 REAL(wp) :: lambda_r !: 1414 REAL(wp) :: mu_r !: 1415 REAL(wp) :: mu_r_2 !: 1416 REAL(wp) :: mu_r_5d2 !: 1417 REAL(wp) :: nr_0 !: 1418 REAL(wp) :: q_s !: 1419 REAL(wp) :: sat !: 1420 REAL(wp) :: t_l !: 1421 REAL(wp) :: temp !: 1422 REAL(wp) :: xr !: 1423 1424 DO k = nzb_s_inner(j,i)+1, nzt 1425 IF ( qr_1d(k) > eps_sb ) THEN 1426 ! 1427 !-- Actual liquid water temperature: 1428 t_l = t_d_pt(k) * pt_1d(k) 1429 ! 1430 !-- Saturation vapor pressure at t_l: 1431 e_s = 610.78_wp * EXP( 17.269_wp * ( t_l - 273.16_wp ) / & 1432 ( t_l - 35.86_wp ) & 1433 ) 1434 ! 1435 !-- Computation of saturation humidity: 1436 q_s = 0.622_wp * e_s / ( hyp(k) - 0.378_wp * e_s ) 1437 alpha = 0.622_wp * l_d_r * l_d_cp / ( t_l * t_l ) 1438 q_s = q_s * ( 1.0_wp + alpha * q_1d(k) ) / ( 1.0_wp + alpha * q_s ) 1439 ! 1440 !-- Supersaturation: 1441 sat = ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0_wp 1442 ! 1443 !-- Evaporation needs only to be calculated in subsaturated regions 1444 IF ( sat < 0.0_wp ) THEN 1445 ! 1446 !-- Actual temperature: 1447 temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) ) 1448 1449 g_evap = 1.0_wp / ( ( l_v / ( r_v * temp ) - 1.0_wp ) * l_v / & 1450 ( thermal_conductivity_l * temp ) + & 1451 r_v * temp / ( diff_coeff_l * e_s ) & 1452 ) 1453 ! 1454 !-- Mean weight of rain drops 1455 xr = hyrho(k) * qr_1d(k) / nr_1d(k) 1456 ! 1457 !-- Weight averaged diameter of rain drops: 1458 dr = ( xr * dpirho_l )**( 1.0_wp / 3.0_wp ) 1459 ! 1460 !-- Compute ventilation factor and intercept parameter 1461 !-- (Seifert and Beheng, 2006; Seifert, 2008): 1462 IF ( ventilation_effect ) THEN 1463 ! 1464 !-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005; 1465 !-- Stevens and Seifert, 2008): 1466 mu_r = 10.0_wp * ( 1.0_wp + TANH( 1.2E3_wp * ( dr - 1.4E-3_wp ) ) ) 1467 ! 1468 !-- Slope parameter of gamma distribution (Seifert, 2008): 1469 lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * & 1470 ( mu_r + 1.0_wp ) & 1471 )**( 1.0_wp / 3.0_wp ) / dr 1472 1473 mu_r_2 = mu_r + 2.0_wp 1474 mu_r_5d2 = mu_r + 2.5_wp 1475 1476 f_vent = a_vent * gamm( mu_r_2 ) * lambda_r**( -mu_r_2 ) + & 1477 b_vent * schmidt_p_1d3 * & 1478 SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) * & 1479 lambda_r**( -mu_r_5d2 ) * & 1480 ( 1.0_wp - & 1481 0.5_wp * ( b_term / a_term ) * & 1482 ( lambda_r / ( c_term + lambda_r ) & 1483 )**mu_r_5d2 - & 1484 0.125_wp * ( b_term / a_term )**2 * & 1485 ( lambda_r / ( 2.0_wp * c_term + lambda_r ) & 1486 )**mu_r_5d2 - & 1487 0.0625_wp * ( b_term / a_term )**3 * & 1488 ( lambda_r / ( 3.0_wp * c_term + lambda_r ) & 1489 )**mu_r_5d2 - & 1490 0.0390625_wp * ( b_term / a_term )**4 * & 1491 ( lambda_r / ( 4.0_wp * c_term + lambda_r ) & 1492 )**mu_r_5d2 & 1493 ) 1494 1495 nr_0 = nr_1d(k) * lambda_r**( mu_r + 1.0_wp ) / & 1496 gamm( mu_r + 1.0_wp ) 1497 ELSE 1498 f_vent = 1.0_wp 1499 nr_0 = nr_1d(k) * dr 1500 ENDIF 1501 ! 1502 !-- Evaporation rate of rain water content (Seifert and Beheng, 2006): 1503 evap = 2.0_wp * pi * nr_0 * g_evap * f_vent * sat / hyrho(k) 1504 evap = MAX( evap, -qr_1d(k) / dt_micro ) 1505 evap_nr = MAX( c_evap * evap / xr * hyrho(k), & 1506 -nr_1d(k) / dt_micro ) 1507 1508 qr_1d(k) = qr_1d(k) + evap * dt_micro 1509 nr_1d(k) = nr_1d(k) + evap_nr * dt_micro 1510 1511 ENDIF 1512 ENDIF 1513 1514 ENDDO 1515 1516 END SUBROUTINE evaporation_rain_ij 1517 1518 1519 SUBROUTINE sedimentation_cloud_ij( i, j ) 1520 1521 USE arrays_3d, & 1522 ONLY: ddzu, dzu, nc_1d, pt_1d, q_1d, qc_1d 1523 1524 USE cloud_parameters, & 1525 ONLY: eps_sb, hyrho, l_d_cp, pt_d_t, sed_qc_const 1526 1527 USE constants, & 1528 ONLY: pi 1529 1530 USE control_parameters, & 1531 ONLY: dt_do2d_xy, dt_micro, intermediate_timestep_count 1532 1533 USE indices, & 1534 ONLY: nzb, nzb_s_inner, nzt 1535 1536 USE kinds 1537 1538 IMPLICIT NONE 1539 1540 INTEGER(iwp) :: i !: 1541 INTEGER(iwp) :: j !: 1542 INTEGER(iwp) :: k !: 1543 1544 REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qc !: 1545 1546 ! 1547 !-- Sedimentation of cloud droplets (Ackermann et al., 2009, MWR): 1548 sed_qc(nzt+1) = 0.0_wp 1549 1550 DO k = nzt, nzb_s_inner(j,i)+1, -1 1551 IF ( qc_1d(k) > eps_sb ) THEN 1552 sed_qc(k) = sed_qc_const * nc_1d(k)**( -2.0_wp / 3.0_wp ) * & 1553 ( qc_1d(k) * hyrho(k) )**( 5.0_wp / 3.0_wp ) 1554 ELSE 1555 sed_qc(k) = 0.0_wp 1556 ENDIF 1557 1558 sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) / & 1559 dt_micro + sed_qc(k+1) & 1560 ) 1561 1562 q_1d(k) = q_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / & 1563 hyrho(k) * dt_micro 1564 qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / & 1565 hyrho(k) * dt_micro 1566 pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / & 1567 hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro 1568 1569 ENDDO 1570 1571 END SUBROUTINE sedimentation_cloud_ij 1572 1573 1574 SUBROUTINE sedimentation_rain_ij( i, j ) 1575 1576 USE arrays_3d, & 1577 ONLY: ddzu, dzu, nr_1d, pt_1d, q_1d, qr_1d 1578 1579 USE cloud_parameters, & 1580 ONLY: a_term, b_term, c_term, cof, dpirho_l, eps_sb, hyrho, & 1581 limiter_sedimentation, l_d_cp, precipitation_amount, prr, & 1582 pt_d_t, stp 1583 1584 USE control_parameters, & 1585 ONLY: call_microphysics_at_all_substeps, dt_do2d_xy, dt_micro, & 1586 dt_3d, intermediate_timestep_count, & 1587 intermediate_timestep_count_max, & 1588 precipitation_amount_interval, time_do2d_xy 1589 1590 USE indices, & 1591 ONLY: nzb, nzb_s_inner, nzt 1592 1593 USE kinds 1594 1595 USE statistics, & 1596 ONLY: weight_substep 1597 1598 IMPLICIT NONE 1599 1600 INTEGER(iwp) :: i !: 1601 INTEGER(iwp) :: j !: 1602 INTEGER(iwp) :: k !: 1603 INTEGER(iwp) :: k_run !: 1604 1605 REAL(wp) :: c_run !: 1606 REAL(wp) :: d_max !: 1607 REAL(wp) :: d_mean !: 1608 REAL(wp) :: d_min !: 1609 REAL(wp) :: dr !: 1610 REAL(wp) :: dt_sedi !: 1611 REAL(wp) :: flux !: 1612 REAL(wp) :: lambda_r !: 1613 REAL(wp) :: mu_r !: 1614 REAL(wp) :: z_run !: 1615 1616 REAL(wp), DIMENSION(nzb:nzt+1) :: c_nr !: 1617 REAL(wp), DIMENSION(nzb:nzt+1) :: c_qr !: 1618 REAL(wp), DIMENSION(nzb:nzt+1) :: d_nr !: 1619 REAL(wp), DIMENSION(nzb:nzt+1) :: d_qr !: 1620 REAL(wp), DIMENSION(nzb:nzt+1) :: nr_slope !: 1621 REAL(wp), DIMENSION(nzb:nzt+1) :: qr_slope !: 1622 REAL(wp), DIMENSION(nzb:nzt+1) :: sed_nr !: 1623 REAL(wp), DIMENSION(nzb:nzt+1) :: sed_qr !: 1624 REAL(wp), DIMENSION(nzb:nzt+1) :: w_nr !: 1625 REAL(wp), DIMENSION(nzb:nzt+1) :: w_qr !: 1626 1627 1628 ! 1629 !-- Computation of sedimentation flux. Implementation according to Stevens 1630 !-- and Seifert (2008). Code is based on UCLA-LES. 948 1631 IF ( intermediate_timestep_count == 1 ) prr(:,j,i) = 0.0_wp 949 1632 ! … … 960 1643 ! 961 1644 !-- Slope parameter of gamma distribution (Seifert, 2008): 962 lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * & 963 ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr 964 965 w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, a_term - b_term * ( 1.0_wp + & 966 c_term / lambda_r )**( -1.0_wp * ( mu_r + 1.0_wp ) ) ) ) 967 w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, a_term - b_term * ( 1.0_wp + & 968 c_term / lambda_r )**( -1.0_wp * ( mu_r + 4.0_wp ) ) ) ) 1645 lambda_r = ( ( mu_r + 3.0_wp ) * ( mu_r + 2.0_wp ) * & 1646 ( mu_r + 1.0_wp ) )**( 1.0_wp / 3.0_wp ) / dr 1647 1648 w_nr(k) = MAX( 0.1_wp, MIN( 20.0_wp, & 1649 a_term - b_term * ( 1.0_wp + & 1650 c_term / lambda_r )**( -1.0_wp * & 1651 ( mu_r + 1.0_wp ) ) & 1652 ) & 1653 ) 1654 w_qr(k) = MAX( 0.1_wp, MIN( 20.0_wp, & 1655 a_term - b_term * ( 1.0_wp + & 1656 c_term / lambda_r )**( -1.0_wp * & 1657 ( mu_r + 4.0_wp ) ) & 1658 ) & 1659 ) 969 1660 ELSE 970 1661 w_nr(k) = 0.0_wp … … 981 1672 !-- Compute Courant number 982 1673 DO k = nzb_s_inner(j,i)+1, nzt 983 c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) * &1674 c_nr(k) = 0.25_wp * ( w_nr(k-1) + 2.0_wp * w_nr(k) + w_nr(k+1) ) * & 984 1675 dt_micro * ddzu(k) 985 c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) * &1676 c_qr(k) = 0.25_wp * ( w_qr(k-1) + 2.0_wp * w_qr(k) + w_qr(k+1) ) * & 986 1677 dt_micro * ddzu(k) 987 1678 ENDDO … … 995 1686 d_max = MAX( qr_1d(k+1), qr_1d(k), qr_1d(k-1) ) - qr_1d(k) 996 1687 997 qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, 2.0_wp * d_max, & 998 ABS( d_mean ) ) 1688 qr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, & 1689 2.0_wp * d_max, & 1690 ABS( d_mean ) ) 999 1691 1000 1692 d_mean = 0.5_wp * ( nr_1d(k+1) + nr_1d(k-1) ) … … 1002 1694 d_max = MAX( nr_1d(k+1), nr_1d(k), nr_1d(k-1) ) - nr_1d(k) 1003 1695 1004 nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, 2.0_wp * d_max, & 1005 ABS( d_mean ) ) 1696 nr_slope(k) = SIGN(1.0_wp, d_mean) * MIN ( 2.0_wp * d_min, & 1697 2.0_wp * d_max, & 1698 ABS( d_mean ) ) 1006 1699 ENDDO 1007 1700 … … 1026 1719 c_run = MIN( 1.0_wp, c_nr(k) ) 1027 1720 DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt ) 1028 flux = flux + hyrho(k_run) * &1029 ( nr_1d(k_run) + nr_slope(k_run) * ( 1.0_wp - c_run ) * 1721 flux = flux + hyrho(k_run) * & 1722 ( nr_1d(k_run) + nr_slope(k_run) * ( 1.0_wp - c_run ) * & 1030 1723 0.5_wp ) * c_run * dzu(k_run) 1031 1724 z_run = z_run + dzu(k_run) … … 1036 1729 !-- It is not allowed to sediment more rain drop number density than 1037 1730 !-- available 1038 flux = MIN( flux, &1731 flux = MIN( flux, & 1039 1732 hyrho(k) * dzu(k+1) * nr_1d(k) + sed_nr(k+1) * dt_micro ) 1040 1733 1041 1734 sed_nr(k) = flux / dt_micro 1042 nr_1d(k) = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) / &1043 hyrho(k) * dt_micro1735 nr_1d(k) = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) / & 1736 hyrho(k) * dt_micro 1044 1737 ! 1045 1738 !-- Sum up all rain water content which contributes to the flux … … 1050 1743 c_run = MIN( 1.0_wp, c_qr(k) ) 1051 1744 1052 DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt -1)1053 1054 flux = flux + hyrho(k_run) * &1055 ( qr_1d(k_run) + qr_slope(k_run) * ( 1.0_wp - c_run ) * 1745 DO WHILE ( c_run > 0.0_wp .AND. k_run <= nzt ) 1746 1747 flux = flux + hyrho(k_run) * & 1748 ( qr_1d(k_run) + qr_slope(k_run) * ( 1.0_wp - c_run ) * & 1056 1749 0.5_wp ) * c_run * dzu(k_run) 1057 1750 z_run = z_run + dzu(k_run) … … 1062 1755 ! 1063 1756 !-- It is not allowed to sediment more rain water content than available 1064 flux = MIN( flux, &1757 flux = MIN( flux, & 1065 1758 hyrho(k) * dzu(k) * qr_1d(k) + sed_qr(k+1) * dt_micro ) 1066 1759 1067 1760 sed_qr(k) = flux / dt_micro 1068 1761 1069 qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &1762 qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / & 1070 1763 hyrho(k) * dt_micro 1071 q_1d(k) = q_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &1764 q_1d(k) = q_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / & 1072 1765 hyrho(k) * dt_micro 1073 pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / &1766 pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / & 1074 1767 hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro 1075 1768 ! 1076 1769 !-- Compute the rain rate 1077 prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * & 1078 weight_substep(intermediate_timestep_count) 1770 IF ( call_microphysics_at_all_substeps ) THEN 1771 prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * & 1772 weight_substep(intermediate_timestep_count) 1773 ELSE 1774 prr(k,j,i) = sed_qr(k) / hyrho(k) 1775 ENDIF 1776 1079 1777 ENDDO 1080 1778 1081 1779 ! 1082 1780 !-- Precipitation amount 1083 IF ( intermediate_timestep_count == intermediate_timestep_count_max &1084 .AND. ( dt_do2d_xy - time_do2d_xy ) < &1085 precipitation_amount_interval ) THEN1086 1087 precipitation_amount(j,i) = precipitation_amount(j,i) + &1088 prr(nzb_s_inner(j,i)+1,j,i) * &1781 IF ( intermediate_timestep_count == intermediate_timestep_count_max & 1782 .AND. ( dt_do2d_xy - time_do2d_xy ) < & 1783 precipitation_amount_interval ) THEN 1784 1785 precipitation_amount(j,i) = precipitation_amount(j,i) + & 1786 prr(nzb_s_inner(j,i)+1,j,i) * & 1089 1787 hyrho(nzb_s_inner(j,i)+1) * dt_3d 1090 1788 ENDIF … … 1092 1790 END SUBROUTINE sedimentation_rain_ij 1093 1791 1094 1792 !------------------------------------------------------------------------------! 1793 ! Call for all optimizations 1794 !------------------------------------------------------------------------------! 1095 1795 ! 1096 1796 !-- This function computes the gamma function (Press et al., 1992). -
palm/trunk/SOURCE/modules.f90
r1360 r1361 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! tend_* removed 23 ! call_microphysics_at_all_substeps added 24 ! default of drizzle set to true 23 25 ! 24 26 ! Former revisions: … … 287 289 diss_l_v, diss_l_w, flux_l_e, flux_l_nr, flux_l_pt, flux_l_q, & 288 290 flux_l_qr, flux_l_sa, flux_l_u, flux_l_v, flux_l_w, kh, km, lad_s, & 289 lad_u, lad_v, lad_w, lai, l_wall, p_loc, sec, sls, tend, t end_pt,&290 tend_nr, tend_q, tend_qr, tric, u_m_l, u_m_n, u_m_r, u_m_s, v_m_l,&291 v_m_n, v_m_r, v_m_s, w_m_l,w_m_n, w_m_r, w_m_s291 lad_u, lad_v, lad_w, lai, l_wall, p_loc, sec, sls, tend, tric, & 292 u_m_l, u_m_n, u_m_r, u_m_s, v_m_l, v_m_n, v_m_r, v_m_s, w_m_l, & 293 w_m_n, w_m_r, w_m_s 292 294 293 295 #if defined( __nopointer ) … … 368 370 b_vent = 0.308_wp, & !: coef. for ventilation effect 369 371 beta_cc = 3.09E-4_wp, & !: coef. in turb. parametrization (cm-2 s3) 370 bfactor, &372 bfactor, & !: solution effect on diffusional growth 371 373 c_1 = 4.82E-6_wp, & !: coef. in turb. parametrization (m) 372 374 c_2 = 4.8E-6_wp, & !: coef. in turb. parametrization (m) … … 384 386 cp = 1005.0_wp, & !: heat capacity of dry air (J kg-1 K-1) 385 387 diff_coeff_l = 0.23E-4_wp, & !: diffusivity of water vapor (m2 s-1) 386 effective_coll_efficiency, & !: 388 effective_coll_efficiency, & !: effective collision efficiency 387 389 eps_ros = 1.0E-4_wp, & !: accuracy of Rosenbrock method 388 390 eps_sb = 1.0E-20_wp, & !: threshold in two-moments scheme … … 408 410 schmidt = 0.71_wp, & !: Schmidt number 409 411 schmidt_p_1d3, & !: schmidt**( 1.0 / 3.0 ) 412 sed_qc_const, & !: const. for sedimentation of cloud water 410 413 sigma_gc = 1.3_wp, & !: log-normal geometric standard deviation 411 414 stp = 2.5066282746310005_wp, & !: parameter in gamma function … … 592 595 bc_lr_raddir = .FALSE., bc_ns_cyc = .TRUE., & 593 596 bc_ns_dirrad = .FALSE., bc_ns_raddir = .FALSE.,& 597 call_microphysics_at_all_substeps = .FALSE., & 594 598 call_psolver_at_all_substeps = .TRUE., & 595 599 cloud_droplets = .FALSE., cloud_physics = .FALSE., & … … 604 608 do_sum = .FALSE., & 605 609 dp_external = .FALSE., dp_smooth = .FALSE., & 606 drizzle = . FALSE., dt_fixed = .FALSE., &610 drizzle = .TRUE., dt_fixed = .FALSE., & 607 611 dt_3d_reached, dt_3d_reached_l, exchange_mg = .FALSE., & 608 612 first_call_lpm = .TRUE., & -
palm/trunk/SOURCE/parin.f90
r1360 r1361 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! +call_microphysics_at_all_substeps 23 23 ! 24 24 ! Former revisions: … … 148 148 bottom_salinityflux, building_height, building_length_x, & 149 149 building_length_y, building_wall_left, building_wall_south, & 150 call_psolver_at_all_substeps, canopy_mode, canyon_height, & 150 call_microphysics_at_all_substeps, call_psolver_at_all_substeps,& 151 canopy_mode, canyon_height, & 151 152 canyon_width_x, canyon_width_y, canyon_wall_left, & 152 153 canyon_wall_south, cfl_factor, & … … 248 249 bottom_salinityflux, building_height, building_length_x, & 249 250 building_length_y, building_wall_left, building_wall_south, & 250 call_psolver_at_all_substeps, canopy_mode, canyon_height, & 251 call_psolver_at_all_substeps, call_microphysics_at_all_substeps, & 252 canopy_mode, canyon_height, & 251 253 canyon_width_x, canyon_width_y, canyon_wall_left, & 252 254 canyon_wall_south, c_sedimentation, cfl_factor, cloud_droplets, & -
palm/trunk/SOURCE/prandtl_fluxes.f90
r1341 r1361 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Bugfix: calculation of turbulent fluxes of rain water content (qrsws) and rain 23 ! drop concentration (nrsws) added 23 24 ! 24 25 ! Former revisions: … … 64 65 65 66 USE arrays_3d, & 66 ONLY: e, pt, q, qs, qsws, rif, shf, ts, u, us, usws, v, vpt, vsws,&67 zu, zw, z0, z0h67 ONLY: e, nr, nrs, nrsws, pt, q, qr, qrs, qrsws, qs, qsws, rif, shf, & 68 ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h 68 69 69 70 USE control_parameters, & 70 ONLY: constant_heatflux, constant_waterflux, coupling_mode, g, & 71 humidity, ibc_e_b, kappa, large_scale_forcing, lsf_surf, & 72 passive_scalar, pt_surface, q_surface, rif_max, rif_min, & 73 run_coupled, surface_pressure 71 ONLY: cloud_physics, constant_heatflux, constant_waterflux, & 72 coupling_mode, g, humidity, ibc_e_b, icloud_scheme, kappa, & 73 large_scale_forcing, lsf_surf, passive_scalar, precipitation, & 74 pt_surface, q_surface, rif_max, rif_min, run_coupled, & 75 surface_pressure 74 76 75 77 USE indices, & … … 96 98 ! 97 99 !-- Data information for accelerators 98 !$acc data present( e, nzb_u_inner, nzb_v_inner, nzb_s_inner, pt, q, qs ) & 99 !$acc present( qsws, rif, shf, ts, u, us, usws, v, vpt, vsws, zu, zw, z0, z0h ) 100 !$acc data present( e, nrsws, nzb_u_inner, nzb_v_inner, nzb_s_inner, pt ) & 101 !$acc present( q, qs, qsws, qrsws, rif, shf, ts, u, us, usws, v ) & 102 !$acc present( vpt, vsws, zu, zw, z0, z0h ) 100 103 ! 101 104 !-- Compute theta* … … 383 386 ENDDO 384 387 ENDIF 388 389 IF ( cloud_physics .AND. icloud_scheme == 0 & 390 .AND. precipitation ) THEN 391 392 !$OMP PARALLEL DO PRIVATE( a, b, k, z_p ) 393 !$acc kernels loop independent 394 DO i = nxlg, nxrg 395 !$acc loop independent 396 DO j = nysg, nyng 397 398 k = nzb_s_inner(j,i) 399 z_p = zu(k+1) - zw(k) 400 401 IF ( rif(j,i) >= 0.0 ) THEN 402 ! 403 !-- Stable stratification 404 qrs(j,i) = kappa * ( qr(k+1,j,i) - qr(k,j,i) ) / ( & 405 LOG( z_p / z0h(j,i) ) + & 406 5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p ) 407 nrs(j,i) = kappa * ( nr(k+1,j,i) - nr(k,j,i) ) / ( & 408 LOG( z_p / z0h(j,i) ) + & 409 5.0 * rif(j,i) * ( z_p - z0h(j,i) ) / z_p ) 410 411 ELSE 412 ! 413 !-- Unstable stratification 414 a = SQRT( 1.0 - 16.0 * rif(j,i) ) 415 b = SQRT( 1.0 - 16.0 * rif(j,i) * z0h(j,i) / z_p ) 416 417 qrs(j,i) = kappa * ( qr(k+1,j,i) - qr(k,j,i) ) / ( & 418 LOG( z_p / z0h(j,i) ) - & 419 2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) ) 420 nrs(j,i) = kappa * ( nr(k+1,j,i) - nr(k,j,i) ) / ( & 421 LOG( z_p / z0h(j,i) ) - & 422 2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) ) 423 424 ENDIF 425 426 ENDDO 427 ENDDO 428 429 ENDIF 430 385 431 ENDIF 386 432 … … 396 442 CALL exchange_horiz_2d( qsws ) 397 443 !$acc update device( qsws ) 444 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. & 445 precipitation ) THEN 446 !$acc update host( qrsws, nrsws ) 447 CALL exchange_horiz_2d( qrsws ) 448 CALL exchange_horiz_2d( nrsws ) 449 !$acc update device( qrsws, nrsws ) 450 ENDIF 398 451 ENDIF 399 452 … … 425 478 426 479 ! 480 !-- Compute (turbulent) fluxes of rain water content and rain drop concentartion 481 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. precipitation ) THEN 482 !$OMP PARALLEL DO 483 !$acc kernels loop independent 484 DO i = nxlg, nxrg 485 !$acc loop independent 486 DO j = nysg, nyng 487 qrsws(j,i) = -qrs(j,i) * us(j,i) 488 nrsws(j,i) = -nrs(j,i) * us(j,i) 489 ENDDO 490 ENDDO 491 ENDIF 492 493 ! 427 494 !-- Bottom boundary condition for the TKE 428 495 IF ( ibc_e_b == 2 ) THEN -
palm/trunk/SOURCE/prognostic_equations.f90
r1354 r1361 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! Two-moment microphysics moved to the start of prognostic equations. This makes 23 ! the 3d arrays for tend_q, tend_qr, tend_pt and tend_pt redundant. 24 ! Additionally, it is allowed to call the microphysics just once during the time 25 ! step (not at each sub-time step). 26 ! 27 ! Two-moment cloud physics added for vector and accelerator optimization. 23 28 ! 24 29 ! Former revisions: … … 126 131 nrswst, pt, ptdf_x, ptdf_y, pt_init, pt_p, prho, q, q_init, & 127 132 q_p, qsws, qswst, qr, qr_p, qrsws, qrswst, rdf, rdf_sc, rho, & 128 sa, sa_init, sa_p, saswsb, saswst, shf, tend, tend_nr, & 129 tend_pt, tend_q, tend_qr, te_m, tnr_m, tpt_m, tq_m, tqr_m, & 130 tsa_m, tswst, tu_m, tv_m, tw_m, u, ug, u_p, v, vg, vpt, v_p, & 131 w, w_p 133 sa, sa_init, sa_p, saswsb, saswst, shf, tend, te_m, tnr_m, & 134 tpt_m, tq_m, tqr_m, tsa_m, tswst, tu_m, tv_m, tw_m, u, ug, u_p, & 135 v, vg, vpt, v_p, w, w_p 132 136 133 137 USE control_parameters, & 134 ONLY: cloud_physics, constant_diffusion, cthf, dp_external, & 138 ONLY: call_microphysics_at_all_substeps, cloud_physics, & 139 constant_diffusion, cthf, dp_external, & 135 140 dp_level_ind_b, dp_smooth_factor, dpdxy, dt_3d, humidity, & 136 141 icloud_scheme, inflow_l, intermediate_timestep_count, & … … 301 306 302 307 DO j = nys, nyn 308 ! 309 !-- If required, calculate cloud microphysical impacts (two-moment scheme) 310 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. & 311 ( intermediate_timestep_count == 1 .OR. & 312 call_microphysics_at_all_substeps ) & 313 ) THEN 314 CALL microphysics_control( i, j ) 315 ENDIF 303 316 ! 304 317 !-- Tendency terms for u-velocity component … … 481 494 ENDIF 482 495 ENDIF 483 ! 484 !-- If required, calculate tendencies for total water content, liquid water 485 !-- potential temperature, rain water content and rain drop concentration 486 IF ( cloud_physics .AND. icloud_scheme == 0 ) CALL microphysics_control( i, j ) 496 487 497 ! 488 498 !-- If required, compute prognostic equation for potential temperature … … 511 521 512 522 ! 513 !-- Using microphysical tendencies (latent heat) 514 IF ( cloud_physics ) THEN 515 IF ( icloud_scheme == 0 ) THEN 516 tend(:,j,i) = tend(:,j,i) + tend_pt(:,j,i) 517 ELSEIF ( icloud_scheme == 1 .AND. precipitation ) THEN 518 CALL impact_of_latent_heat( i, j ) 519 ENDIF 523 !-- If required compute impact of latent heat due to precipitation 524 IF ( cloud_physics .AND. icloud_scheme == 1 .AND. & 525 precipitation ) THEN 526 CALL impact_of_latent_heat( i, j ) 520 527 ENDIF 521 528 … … 641 648 642 649 ! 643 !-- Using microphysical tendencies 644 IF ( cloud_physics ) THEN 645 IF ( icloud_scheme == 0 ) THEN 646 tend(:,j,i) = tend(:,j,i) + tend_q(:,j,i) 647 ELSEIF ( icloud_scheme == 1 .AND. precipitation ) THEN 648 CALL calc_precipitation( i, j ) 649 ENDIF 650 !-- If required compute decrease of total water content due to 651 !-- precipitation 652 IF ( cloud_physics .AND. icloud_scheme == 1 .AND. & 653 precipitation ) THEN 654 CALL calc_precipitation( i, j ) 650 655 ENDIF 651 656 ! … … 712 717 ENDIF 713 718 CALL diffusion_s( i, j, qr, qrsws, qrswst, wall_qrflux ) 714 !715 !-- Using microphysical tendencies (autoconversion, accretion,716 !-- evaporation; if required: sedimentation)717 tend(:,j,i) = tend(:,j,i) + tend_qr(:,j,i)718 719 719 720 CALL user_actions( i, j, 'qr-tendency' ) … … 757 758 ENDIF 758 759 CALL diffusion_s( i, j, nr, nrsws, nrswst, wall_nrflux ) 759 !760 !-- Using microphysical tendencies (autoconversion, accretion,761 !-- selfcollection, breakup, evaporation;762 !-- if required: sedimentation)763 tend(:,j,i) = tend(:,j,i) + tend_nr(:,j,i)764 760 765 761 CALL user_actions( i, j, 'nr-tendency' ) … … 883 879 884 880 ! 881 !-- If required, calculate cloud microphysical impacts (two-moment scheme) 882 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. & 883 ( intermediate_timestep_count == 1 .OR. & 884 call_microphysics_at_all_substeps ) & 885 ) THEN 886 CALL cpu_log( log_point(51), 'microphysics', 'start' ) 887 CALL microphysics_control 888 CALL cpu_log( log_point(51), 'microphysics', 'stop' ) 889 ENDIF 890 891 ! 885 892 !-- u-velocity component 886 893 CALL cpu_log( log_point(5), 'u-equation', 'start' ) … … 1156 1163 ! 1157 1164 !-- If required compute impact of latent heat due to precipitation 1158 IF ( precipitation ) THEN1165 IF ( cloud_physics .AND. icloud_scheme == 1 .AND. precipitation ) THEN 1159 1166 CALL impact_of_latent_heat 1160 1167 ENDIF … … 1348 1355 !-- If required compute decrease of total water content due to 1349 1356 !-- precipitation 1350 IF ( precipitation ) THEN1357 IF ( cloud_physics .AND. icloud_scheme == 1 .AND. precipitation ) THEN 1351 1358 CALL calc_precipitation 1352 1359 ENDIF … … 1406 1413 1407 1414 CALL cpu_log( log_point(29), 'q/s-equation', 'stop' ) 1415 1416 ! 1417 !-- If required, calculate prognostic equations for rain water content 1418 !-- and rain drop concentration 1419 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. precipitation ) THEN 1420 1421 CALL cpu_log( log_point(52), 'qr-equation', 'start' ) 1422 1423 ! 1424 !-- Calculate prognostic equation for rain water content 1425 sbt = tsc(2) 1426 IF ( scalar_advec == 'bc-scheme' ) THEN 1427 1428 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 1429 ! 1430 !-- Bott-Chlond scheme always uses Euler time step. Thus: 1431 sbt = 1.0_wp 1432 ENDIF 1433 tend = 0.0_wp 1434 CALL advec_s_bc( qr, 'qr' ) 1435 1436 ENDIF 1437 1438 ! 1439 !-- qr-tendency terms with no communication 1440 IF ( scalar_advec /= 'bc-scheme' ) THEN 1441 tend = 0.0_wp 1442 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1443 IF ( ws_scheme_sca ) THEN 1444 CALL advec_s_ws( qr, 'qr' ) 1445 ELSE 1446 CALL advec_s_pw( qr ) 1447 ENDIF 1448 ELSE 1449 CALL advec_s_up( qr ) 1450 ENDIF 1451 ENDIF 1452 1453 CALL diffusion_s( qr, qrsws, qrswst, wall_qrflux ) 1454 1455 CALL user_actions( 'qr-tendency' ) 1456 1457 ! 1458 !-- Prognostic equation for rain water content 1459 DO i = nxl, nxr 1460 DO j = nys, nyn 1461 DO k = nzb_s_inner(j,i)+1, nzt 1462 qr_p(k,j,i) = qr(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + & 1463 tsc(3) * tqr_m(k,j,i) ) & 1464 - tsc(5) * rdf_sc(k) * qr(k,j,i) 1465 IF ( qr_p(k,j,i) < 0.0_wp ) qr_p(k,j,i) = 0.0_wp 1466 ENDDO 1467 ENDDO 1468 ENDDO 1469 1470 ! 1471 !-- Calculate tendencies for the next Runge-Kutta step 1472 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1473 IF ( intermediate_timestep_count == 1 ) THEN 1474 DO i = nxl, nxr 1475 DO j = nys, nyn 1476 DO k = nzb_s_inner(j,i)+1, nzt 1477 tqr_m(k,j,i) = tend(k,j,i) 1478 ENDDO 1479 ENDDO 1480 ENDDO 1481 ELSEIF ( intermediate_timestep_count < & 1482 intermediate_timestep_count_max ) THEN 1483 DO i = nxl, nxr 1484 DO j = nys, nyn 1485 DO k = nzb_s_inner(j,i)+1, nzt 1486 tqr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * & 1487 tqr_m(k,j,i) 1488 ENDDO 1489 ENDDO 1490 ENDDO 1491 ENDIF 1492 ENDIF 1493 1494 CALL cpu_log( log_point(52), 'qr-equation', 'stop' ) 1495 CALL cpu_log( log_point(53), 'nr-equation', 'start' ) 1496 1497 ! 1498 !-- Calculate prognostic equation for rain drop concentration 1499 sbt = tsc(2) 1500 IF ( scalar_advec == 'bc-scheme' ) THEN 1501 1502 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 1503 ! 1504 !-- Bott-Chlond scheme always uses Euler time step. Thus: 1505 sbt = 1.0_wp 1506 ENDIF 1507 tend = 0.0_wp 1508 CALL advec_s_bc( nr, 'nr' ) 1509 1510 ENDIF 1511 1512 ! 1513 !-- nr-tendency terms with no communication 1514 IF ( scalar_advec /= 'bc-scheme' ) THEN 1515 tend = 0.0_wp 1516 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1517 IF ( ws_scheme_sca ) THEN 1518 CALL advec_s_ws( nr, 'nr' ) 1519 ELSE 1520 CALL advec_s_pw( nr ) 1521 ENDIF 1522 ELSE 1523 CALL advec_s_up( nr ) 1524 ENDIF 1525 ENDIF 1526 1527 CALL diffusion_s( nr, nrsws, nrswst, wall_nrflux ) 1528 1529 CALL user_actions( 'nr-tendency' ) 1530 1531 ! 1532 !-- Prognostic equation for rain drop concentration 1533 DO i = nxl, nxr 1534 DO j = nys, nyn 1535 DO k = nzb_s_inner(j,i)+1, nzt 1536 nr_p(k,j,i) = nr(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + & 1537 tsc(3) * tnr_m(k,j,i) ) & 1538 - tsc(5) * rdf_sc(k) * nr(k,j,i) 1539 IF ( nr_p(k,j,i) < 0.0_wp ) nr_p(k,j,i) = 0.0_wp 1540 ENDDO 1541 ENDDO 1542 ENDDO 1543 1544 ! 1545 !-- Calculate tendencies for the next Runge-Kutta step 1546 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1547 IF ( intermediate_timestep_count == 1 ) THEN 1548 DO i = nxl, nxr 1549 DO j = nys, nyn 1550 DO k = nzb_s_inner(j,i)+1, nzt 1551 tnr_m(k,j,i) = tend(k,j,i) 1552 ENDDO 1553 ENDDO 1554 ENDDO 1555 ELSEIF ( intermediate_timestep_count < & 1556 intermediate_timestep_count_max ) THEN 1557 DO i = nxl, nxr 1558 DO j = nys, nyn 1559 DO k = nzb_s_inner(j,i)+1, nzt 1560 tnr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * & 1561 tnr_m(k,j,i) 1562 ENDDO 1563 ENDDO 1564 ENDDO 1565 ENDIF 1566 ENDIF 1567 1568 CALL cpu_log( log_point(53), 'nr-equation', 'stop' ) 1569 1570 ENDIF 1408 1571 1409 1572 ENDIF … … 1510 1673 ENDIF 1511 1674 1512 1513 1675 END SUBROUTINE prognostic_equations_vector 1514 1676 … … 1539 1701 runge_step = 2 1540 1702 ENDIF 1703 ENDIF 1704 1705 ! 1706 !-- If required, calculate cloud microphysical impacts (two-moment scheme) 1707 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. & 1708 ( intermediate_timestep_count == 1 .OR. & 1709 call_microphysics_at_all_substeps ) & 1710 ) THEN 1711 CALL cpu_log( log_point(51), 'microphysics', 'start' ) 1712 CALL microphysics_control 1713 CALL cpu_log( log_point(51), 'microphysics', 'stop' ) 1541 1714 ENDIF 1542 1715 … … 1791 1964 ! 1792 1965 !-- If required compute impact of latent heat due to precipitation 1793 IF ( precipitation ) THEN1966 IF ( cloud_physics .AND. icloud_scheme == 1 .AND. precipitation ) THEN 1794 1967 CALL impact_of_latent_heat 1795 1968 ENDIF … … 1957 2130 !-- If required compute decrease of total water content due to 1958 2131 !-- precipitation 1959 IF ( precipitation ) THEN2132 IF ( cloud_physics .AND. icloud_scheme == 1 .AND. precipitation ) THEN 1960 2133 CALL calc_precipitation 1961 2134 ENDIF … … 1999 2172 2000 2173 CALL cpu_log( log_point(29), 'q/s-equation', 'stop' ) 2174 2175 ! 2176 !-- If required, calculate prognostic equations for rain water content 2177 !-- and rain drop concentration 2178 IF ( cloud_physics .AND. icloud_scheme == 0 .AND. precipitation ) THEN 2179 2180 CALL cpu_log( log_point(52), 'qr-equation', 'start' ) 2181 ! 2182 !-- qr-tendency terms with communication 2183 sbt = tsc(2) 2184 IF ( scalar_advec == 'bc-scheme' ) THEN 2185 2186 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 2187 ! 2188 !-- Bott-Chlond scheme always uses Euler time step. Thus: 2189 sbt = 1.0_wp 2190 ENDIF 2191 tend = 0.0_wp 2192 CALL advec_s_bc( qr, 'qr' ) 2193 2194 ENDIF 2195 2196 ! 2197 !-- qr-tendency terms with no communication 2198 IF ( scalar_advec /= 'bc-scheme' ) THEN 2199 tend = 0.0_wp 2200 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2201 IF ( ws_scheme_sca ) THEN 2202 CALL advec_s_ws( qr, 'qr' ) 2203 ELSE 2204 CALL advec_s_pw( qr ) 2205 ENDIF 2206 ELSE 2207 CALL advec_s_up( qr ) 2208 ENDIF 2209 ENDIF 2210 2211 CALL diffusion_s( qr, qrsws, qrswst, wall_qrflux ) 2212 2213 CALL user_actions( 'qr-tendency' ) 2214 2215 ! 2216 !-- Prognostic equation for rain water content 2217 DO i = i_left, i_right 2218 DO j = j_south, j_north 2219 DO k = nzb_s_inner(j,i)+1, nzt 2220 qr_p(k,j,i) = qr(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + & 2221 tsc(3) * tqr_m(k,j,i) ) & 2222 - tsc(5) * rdf_sc(k) * qr(k,j,i) 2223 IF ( qr_p(k,j,i) < 0.0_wp ) qr_p(k,j,i) = 0.0_wp 2224 ! 2225 !-- Tendencies for the next Runge-Kutta step 2226 IF ( runge_step == 1 ) THEN 2227 tqr_m(k,j,i) = tend(k,j,i) 2228 ELSEIF ( runge_step == 2 ) THEN 2229 tqr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * & 2230 tqr_m(k,j,i) 2231 ENDIF 2232 ENDDO 2233 ENDDO 2234 ENDDO 2235 2236 CALL cpu_log( log_point(52), 'qr-equation', 'stop' ) 2237 CALL cpu_log( log_point(53), 'nr-equation', 'start' ) 2238 2239 ! 2240 !-- nr-tendency terms with communication 2241 sbt = tsc(2) 2242 IF ( scalar_advec == 'bc-scheme' ) THEN 2243 2244 IF ( timestep_scheme(1:5) /= 'runge' ) THEN 2245 ! 2246 !-- Bott-Chlond scheme always uses Euler time step. Thus: 2247 sbt = 1.0_wp 2248 ENDIF 2249 tend = 0.0_wp 2250 CALL advec_s_bc( nr, 'nr' ) 2251 2252 ENDIF 2253 2254 ! 2255 !-- nr-tendency terms with no communication 2256 IF ( scalar_advec /= 'bc-scheme' ) THEN 2257 tend = 0.0_wp 2258 IF ( timestep_scheme(1:5) == 'runge' ) THEN 2259 IF ( ws_scheme_sca ) THEN 2260 CALL advec_s_ws( nr, 'nr' ) 2261 ELSE 2262 CALL advec_s_pw( nr ) 2263 ENDIF 2264 ELSE 2265 CALL advec_s_up( nr ) 2266 ENDIF 2267 ENDIF 2268 2269 CALL diffusion_s( nr, nrsws, nrswst, wall_nrflux ) 2270 2271 CALL user_actions( 'nr-tendency' ) 2272 2273 ! 2274 !-- Prognostic equation for rain drop concentration 2275 DO i = i_left, i_right 2276 DO j = j_south, j_north 2277 DO k = nzb_s_inner(j,i)+1, nzt 2278 nr_p(k,j,i) = nr(k,j,i) + dt_3d * ( sbt * tend(k,j,i) + & 2279 tsc(3) * tnr_m(k,j,i) ) & 2280 - tsc(5) * rdf_sc(k) * nr(k,j,i) 2281 IF ( nr_p(k,j,i) < 0.0_wp ) nr_p(k,j,i) = 0.0_wp 2282 ! 2283 !-- Tendencies for the next Runge-Kutta step 2284 IF ( runge_step == 1 ) THEN 2285 tnr_m(k,j,i) = tend(k,j,i) 2286 ELSEIF ( runge_step == 2 ) THEN 2287 tnr_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * & 2288 tnr_m(k,j,i) 2289 ENDIF 2290 ENDDO 2291 ENDDO 2292 ENDDO 2293 2294 CALL cpu_log( log_point(53), 'nr-equation', 'stop' ) 2295 2296 ENDIF 2001 2297 2002 2298 ENDIF … … 2095 2391 ENDIF 2096 2392 2097 2098 2393 END SUBROUTINE prognostic_equations_acc 2099 2394
Note: See TracChangeset
for help on using the changeset viewer.