Changeset 1115 for palm/trunk/SOURCE/microphysics.f90
- Timestamp:
- Mar 26, 2013 6:16:16 PM (11 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/microphysics.f90 (modified) (40 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/microphysics.f90
r1107 r1115 20 20 ! Current revisions: 21 21 ! ------------------ 22 ! 22 ! microphyical tendencies are calculated in microphysics_control in an optimized 23 ! way; unrealistic values are prevented; bugfix in evaporation; some reformatting 23 24 ! 24 25 ! Former revisions: … … 35 36 ! 1065 2012-11-22 17:42:36Z hoffmann 36 37 ! Sedimentation process implemented according to Stevens and Seifert (2008). 37 ! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens 38 ! Turbulence effects on autoconversion and accretion added (Seifert, Nuijens 38 39 ! and Stevens, 2010). 39 40 ! … … 48 49 49 50 PRIVATE 50 PUBLIC dsd_properties, autoconversion, accretion, selfcollection_breakup, & 51 evaporation_rain, sedimentation_cloud, sedimentation_rain 52 53 INTERFACE dsd_properties 54 MODULE PROCEDURE dsd_properties 55 MODULE PROCEDURE dsd_properties_ij 56 END INTERFACE dsd_properties 51 PUBLIC microphysics_control 52 53 INTERFACE microphysics_control 54 MODULE PROCEDURE microphysics_control 55 MODULE PROCEDURE microphysics_control_ij 56 END INTERFACE microphysics_control 57 58 INTERFACE adjust_cloud 59 MODULE PROCEDURE adjust_cloud 60 MODULE PROCEDURE adjust_cloud_ij 61 END INTERFACE adjust_cloud 57 62 58 63 INTERFACE autoconversion … … 92 97 ! Call for all grid points 93 98 !------------------------------------------------------------------------------! 94 SUBROUTINE dsd_properties95 96 USE arrays_3d 97 USE c loud_parameters98 USE constants99 USE indices99 SUBROUTINE microphysics_control 100 101 USE arrays_3d 102 USE control_parameters 103 USE indices 104 USE statistics 100 105 101 106 IMPLICIT NONE … … 106 111 DO i = nxl, nxr 107 112 DO j = nys, nyn 108 DO k = nzb_ 2d(j,i)+1, nzt113 DO k = nzb_s_inner(j,i)+1, nzt 109 114 110 115 ENDDO … … 112 117 ENDDO 113 118 114 END SUBROUTINE dsd_properties 115 116 117 SUBROUTINE autoconversion 118 119 USE arrays_3d 120 USE cloud_parameters 121 USE constants 119 END SUBROUTINE microphysics_control 120 121 SUBROUTINE adjust_cloud 122 123 USE arrays_3d 124 USE cloud_parameters 122 125 USE indices 123 126 … … 129 132 DO i = nxl, nxr 130 133 DO j = nys, nyn 131 DO k = nzb_ 2d(j,i)+1, nzt134 DO k = nzb_s_inner(j,i)+1, nzt 132 135 133 136 ENDDO … … 135 138 ENDDO 136 139 137 END SUBROUTINE autoconversion 138 139 140 SUBROUTINE accretion 141 142 USE arrays_3d 143 USE cloud_parameters 144 USE constants 140 END SUBROUTINE adjust_cloud 141 142 143 SUBROUTINE autoconversion 144 145 USE arrays_3d 146 USE cloud_parameters 147 USE control_parameters 148 USE grid_variables 145 149 USE indices 146 150 … … 152 156 DO i = nxl, nxr 153 157 DO j = nys, nyn 154 DO k = nzb_ 2d(j,i)+1, nzt158 DO k = nzb_s_inner(j,i)+1, nzt 155 159 156 160 ENDDO … … 158 162 ENDDO 159 163 160 END SUBROUTINE a ccretion161 162 163 SUBROUTINE selfcollection_breakup164 165 USE arrays_3d 166 USE cloud_parameters 167 USE con stants164 END SUBROUTINE autoconversion 165 166 167 SUBROUTINE accretion 168 169 USE arrays_3d 170 USE cloud_parameters 171 USE control_parameters 168 172 USE indices 169 173 … … 175 179 DO i = nxl, nxr 176 180 DO j = nys, nyn 177 DO k = nzb_ 2d(j,i)+1, nzt181 DO k = nzb_s_inner(j,i)+1, nzt 178 182 179 183 ENDDO … … 181 185 ENDDO 182 186 183 END SUBROUTINE selfcollection_breakup184 185 186 SUBROUTINE evaporation_rain187 188 USE arrays_3d 189 USE cloud_parameters 190 USE con stants187 END SUBROUTINE accretion 188 189 190 SUBROUTINE selfcollection_breakup 191 192 USE arrays_3d 193 USE cloud_parameters 194 USE control_parameters 191 195 USE indices 192 196 … … 198 202 DO i = nxl, nxr 199 203 DO j = nys, nyn 200 DO k = nzb_ 2d(j,i)+1, nzt204 DO k = nzb_s_inner(j,i)+1, nzt 201 205 202 206 ENDDO … … 204 208 ENDDO 205 209 206 END SUBROUTINE evaporation_rain207 208 209 SUBROUTINE sedimentation_cloud210 END SUBROUTINE selfcollection_breakup 211 212 213 SUBROUTINE evaporation_rain 210 214 211 215 USE arrays_3d 212 216 USE cloud_parameters 213 217 USE constants 218 USE control_parameters 214 219 USE indices 215 220 … … 221 226 DO i = nxl, nxr 222 227 DO j = nys, nyn 223 DO k = nzb_ 2d(j,i)+1, nzt228 DO k = nzb_s_inner(j,i)+1, nzt 224 229 225 230 ENDDO … … 227 232 ENDDO 228 233 229 END SUBROUTINE sedimentation_cloud230 231 232 SUBROUTINE sedimentation_ rain234 END SUBROUTINE evaporation_rain 235 236 237 SUBROUTINE sedimentation_cloud 233 238 234 239 USE arrays_3d 235 240 USE cloud_parameters 236 241 USE constants 242 USE control_parameters 237 243 USE indices 238 244 … … 244 250 DO i = nxl, nxr 245 251 DO j = nys, nyn 246 DO k = nzb_2d(j,i)+1, nzt 252 DO k = nzb_s_inner(j,i)+1, nzt 253 254 ENDDO 255 ENDDO 256 ENDDO 257 258 END SUBROUTINE sedimentation_cloud 259 260 261 SUBROUTINE sedimentation_rain 262 263 USE arrays_3d 264 USE cloud_parameters 265 USE constants 266 USE control_parameters 267 USE indices 268 USE statistics 269 270 IMPLICIT NONE 271 272 INTEGER :: i, j, k 273 274 275 DO i = nxl, nxr 276 DO j = nys, nyn 277 DO k = nzb_s_inner(j,i)+1, nzt 247 278 248 279 ENDDO … … 256 287 ! Call for grid point i,j 257 288 !------------------------------------------------------------------------------! 258 SUBROUTINE dsd_properties_ij( i, j ) 259 260 USE arrays_3d 261 USE cloud_parameters 262 USE constants 263 USE indices 264 USE control_parameters 265 USE user 266 267 IMPLICIT NONE 268 269 INTEGER :: i, j, k 270 271 DO k = nzb_2d(j,i)+1, nzt 289 290 SUBROUTINE microphysics_control_ij( i, j ) 291 292 USE arrays_3d 293 USE cloud_parameters 294 USE control_parameters 295 USE statistics 296 297 IMPLICIT NONE 298 299 INTEGER :: i, j 300 301 dt_micro = dt_3d * weight_pres(intermediate_timestep_count) 302 ! 303 !-- Adjust unrealistic values 304 IF ( precipitation ) CALL adjust_cloud( i,j ) 305 ! 306 !-- Use 1-d arrays 307 q_1d(:) = q(:,j,i) 308 pt_1d(:) = pt(:,j,i) 309 qc_1d(:) = qc(:,j,i) 310 nc_1d(:) = nc_const 311 IF ( precipitation ) THEN 312 qr_1d(:) = qr(:,j,i) 313 nr_1d(:) = nr(:,j,i) 314 ENDIF 315 ! 316 !-- Compute cloud physics 317 IF ( precipitation ) THEN 318 CALL autoconversion( i,j ) 319 CALL accretion( i,j ) 320 CALL selfcollection_breakup( i,j ) 321 CALL evaporation_rain( i,j ) 322 CALL sedimentation_rain( i,j ) 323 ENDIF 324 325 IF ( drizzle ) CALL sedimentation_cloud( i,j ) 326 ! 327 !-- Derive tendencies 328 tend_q(:,j,i) = ( q_1d(:) - q(:,j,i) ) / dt_micro 329 tend_pt(:,j,i) = ( pt_1d(:) - pt(:,j,i) ) / dt_micro 330 IF ( precipitation ) THEN 331 tend_qr(:,j,i) = ( qr_1d(:) - qr(:,j,i) ) / dt_micro 332 tend_nr(:,j,i) = ( nr_1d(:) - nr(:,j,i) ) / dt_micro 333 ENDIF 334 335 END SUBROUTINE microphysics_control_ij 336 337 SUBROUTINE adjust_cloud_ij( i, j ) 338 339 USE arrays_3d 340 USE cloud_parameters 341 USE indices 342 343 IMPLICIT NONE 344 345 INTEGER :: i, j, k 346 ! 347 !-- Adjust number of raindrops to avoid nonlinear effects in 348 !-- sedimentation and evaporation of rain drops due to too small or 349 !-- too big weights of rain drops (Stevens and Seifert, 2008). 350 !-- The same procedure is applied to cloud droplets if they are determined 351 !-- prognostically. 352 DO k = nzb_s_inner(j,i)+1, nzt 272 353 273 354 IF ( qr(k,j,i) <= eps_sb ) THEN 274 355 qr(k,j,i) = 0.0 356 nr(k,j,i) = 0.0 275 357 ELSE 276 358 ! … … 283 365 nr(k,j,i) = qr(k,j,i) * hyrho(k) / xrmax 284 366 ENDIF 285 xr(k) = hyrho(k) * qr(k,j,i) / nr(k,j,i) 286 ! 287 !-- Weight averaged diameter of rain drops: 288 dr(k) = ( hyrho(k) * qr(k,j,i) / nr(k,j,i) * & 289 dpirho_l )**( 1.0 / 3.0 ) 290 ! 291 !-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005; 292 !-- Stevens and Seifert, 2008): 293 mu_r(k) = 10.0 * ( 1.0 + TANH( 1.2E3 * ( dr(k) - 1.4E-3 ) ) ) 294 ! 295 !-- Slope parameter of gamma distribution (Seifert, 2008): 296 lambda_r(k) = ( ( mu_r(k) + 3.0 ) * ( mu_r(k) + 2.0 ) * & 297 ( mu_r(k) + 1.0 ) )**( 1.0 / 3.0 ) / dr(k) 367 298 368 ENDIF 299 ENDDO 300 301 END SUBROUTINE dsd_properties_ij 369 370 ENDDO 371 372 END SUBROUTINE adjust_cloud_ij 302 373 303 374 … … 306 377 USE arrays_3d 307 378 USE cloud_parameters 308 USE constants 309 USE indices 310 USE control_parameters 311 USE statistics 379 USE control_parameters 312 380 USE grid_variables 313 314 IMPLICIT NONE 315 316 INTEGER :: i, j, k 317 REAL :: k_au, autocon, phi_au, tau_cloud, xc, nu_c, rc, & 318 l_mix, re_lambda, alpha_cc, r_cc, sigma_cc, epsilon 381 USE indices 382 383 IMPLICIT NONE 384 385 INTEGER :: i, j, k 386 REAL :: alpha_cc, autocon, epsilon, k_au, l_mix, nu_c, phi_au, & 387 r_cc, rc, re_lambda, selfcoll, sigma_cc, tau_cloud, xc 319 388 320 389 321 390 k_au = k_cc / ( 20.0 * x0 ) 322 391 323 DO k = nzb_ 2d(j,i)+1, nzt324 325 IF ( q l(k,j,i) > 0.0) THEN392 DO k = nzb_s_inner(j,i)+1, nzt 393 394 IF ( qc_1d(k) > eps_sb ) THEN 326 395 ! 327 396 !-- Intern time scale of coagulation (Seifert and Beheng, 2006): 328 !-- (1.0 - q l(k,j,i) / ( ql(k,j,i) + qr(k,j,i) ))329 tau_cloud = 1.0 - q l(k,j,i) / ( ql(k,j,i) + qr(k,j,i) + 1.0E-20)397 !-- (1.0 - qc(k,j,i) / ( qc(k,j,i) + qr_1d(k) )) 398 tau_cloud = 1.0 - qc_1d(k) / ( qr_1d(k) + qc_1d(k) ) 330 399 ! 331 400 !-- Universal function for autoconversion process … … 335 404 !-- Shape parameter of gamma distribution (Geoffroy et al., 2010): 336 405 !-- (Use constant nu_c = 1.0 instead?) 337 nu_c = 1.0 !MAX( 0.0, 1580.0 * hyrho(k) * q l(k,j,i) - 0.28 )406 nu_c = 1.0 !MAX( 0.0, 1580.0 * hyrho(k) * qc(k,j,i) - 0.28 ) 338 407 ! 339 408 !-- Mean weight of cloud droplets: 340 xc = hyrho(k) * ql(k,j,i) / nc409 xc = hyrho(k) * qc_1d(k) / nc_1d(k) 341 410 ! 342 411 !-- Parameterized turbulence effects on autoconversion (Seifert, … … 371 440 ! 372 441 !-- Autoconversion rate (Seifert and Beheng, 2006): 373 autocon = k_au * ( nu_c + 2.0 ) * ( nu_c + 4.0 ) / & 374 ( nu_c + 1.0 )**2 * ql(k,j,i)**2 * xc**2 * & 375 ( 1.0 + phi_au / ( 1.0 - tau_cloud )**2 ) * & 376 rho_surface 377 autocon = MIN( autocon, ql(k,j,i) / ( dt_3d * & 378 weight_substep(intermediate_timestep_count) ) ) 379 ! 380 !-- Tendencies for q, qr, nr, pt: 381 tend_qr(k,j,i) = tend_qr(k,j,i) + autocon 382 tend_q(k,j,i) = tend_q(k,j,i) - autocon 383 tend_nr(k,j,i) = tend_nr(k,j,i) + autocon / x0 * hyrho(k) 384 tend_pt(k,j,i) = tend_pt(k,j,i) + autocon * l_d_cp * pt_d_t(k) 442 autocon = k_au * ( nu_c + 2.0 ) * ( nu_c + 4.0 ) / & 443 ( nu_c + 1.0 )**2 * qc_1d(k)**2 * xc**2 * & 444 ( 1.0 + phi_au / ( 1.0 - tau_cloud )**2 ) * & 445 rho_surface 446 autocon = MIN( autocon, qc_1d(k) / dt_micro ) 447 448 qr_1d(k) = qr_1d(k) + autocon * dt_micro 449 qc_1d(k) = qc_1d(k) - autocon * dt_micro 450 nr_1d(k) = nr_1d(k) + autocon / x0 * hyrho(k) * dt_micro 385 451 386 452 ENDIF … … 395 461 USE arrays_3d 396 462 USE cloud_parameters 397 USE constants 398 USE indices 399 USE control_parameters 400 USE statistics 401 402 IMPLICIT NONE 403 404 INTEGER :: i, j, k 405 REAL :: accr, phi_ac, tau_cloud, k_cr 406 407 DO k = nzb_2d(j,i)+1, nzt 408 409 IF ( ( ql(k,j,i) > 0.0 ) .AND. ( qr(k,j,i) > eps_sb ) ) THEN 463 USE control_parameters 464 USE indices 465 466 IMPLICIT NONE 467 468 INTEGER :: i, j, k 469 REAL :: accr, k_cr, phi_ac, tau_cloud, xc 470 471 DO k = nzb_s_inner(j,i)+1, nzt 472 IF ( ( qc_1d(k) > eps_sb ) .AND. ( qr_1d(k) > eps_sb ) ) THEN 410 473 ! 411 474 !-- Intern time scale of coagulation (Seifert and Beheng, 2006): 412 tau_cloud = 1.0 - q l(k,j,i) / ( ql(k,j,i) + qr(k,j,i) + 1.0E-20)475 tau_cloud = 1.0 - qc_1d(k) / ( qc_1d(k) + qr_1d(k) ) 413 476 ! 414 477 !-- Universal function for accretion process … … 421 484 !-- convert the dissipation (diss) from m2 s-3 to cm2 s-3. 422 485 IF ( turbulence ) THEN 423 k_cr = k_cr0 * ( 1.0 + 0.05 * &486 k_cr = k_cr0 * ( 1.0 + 0.05 * & 424 487 MIN( 600.0, diss(k,j,i) * 1.0E4 )**0.25 ) 425 488 ELSE … … 428 491 ! 429 492 !-- Accretion rate (Seifert and Beheng, 2006): 430 accr = k_cr * q l(k,j,i) * qr(k,j,i) * phi_ac *&493 accr = k_cr * qc_1d(k) * qr_1d(k) * phi_ac * & 431 494 SQRT( rho_surface * hyrho(k) ) 432 accr = MIN( accr, ql(k,j,i) / ( dt_3d * & 433 weight_substep(intermediate_timestep_count) ) ) 434 ! 435 !-- Tendencies for q, qr, pt: 436 tend_qr(k,j,i) = tend_qr(k,j,i) + accr 437 tend_q(k,j,i) = tend_q(k,j,i) - accr 438 tend_pt(k,j,i) = tend_pt(k,j,i) + accr * l_d_cp * pt_d_t(k) 495 accr = MIN( accr, qc_1d(k) / dt_micro ) 496 497 qr_1d(k) = qr_1d(k) + accr * dt_micro 498 qc_1d(k) = qc_1d(k) - accr * dt_micro 439 499 440 500 ENDIF … … 449 509 USE arrays_3d 450 510 USE cloud_parameters 451 USE constants 452 USE indices 453 USE control_parameters 454 USE statistics 511 USE control_parameters 512 USE indices 455 513 456 514 IMPLICIT NONE 457 515 458 516 INTEGER :: i, j, k 459 REAL :: selfcoll, breakup, phi_br, phi_sc 460 461 462 DO k = nzb_2d(j,i)+1, nzt 463 464 IF ( qr(k,j,i) > eps_sb ) THEN 465 ! 466 !-- Selfcollection rate (Seifert and Beheng, 2006): 467 !-- pirho_l**( 1.0 / 3.0 ) is necessary to convert [lambda_r] = m-1 to 468 !-- kg**( 1.0 / 3.0 ). 469 phi_sc = 1.0 !( 1.0 + kappa_rr / lambda_r(k) * & 470 !pirho_l**( 1.0 / 3.0 ) )**( -9 ) 471 472 selfcoll = k_rr * nr(k,j,i) * qr(k,j,i) * phi_sc * & 517 REAL :: breakup, dr, phi_br, selfcoll 518 519 DO k = nzb_s_inner(j,i)+1, nzt 520 IF ( qr_1d(k) > eps_sb ) THEN 521 ! 522 !-- Selfcollection rate (Seifert and Beheng, 2001): 523 selfcoll = k_rr * nr_1d(k) * qr_1d(k) * & 473 524 SQRT( hyrho(k) * rho_surface ) 474 525 ! 526 !-- Weight averaged diameter of rain drops: 527 dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0 / 3.0 ) 528 ! 475 529 !-- Collisional breakup rate (Seifert, 2008): 476 IF ( dr (k)>= 0.3E-3 ) THEN477 phi_br = k_br * ( dr (k)- 1.1E-3 )530 IF ( dr >= 0.3E-3 ) THEN 531 phi_br = k_br * ( dr - 1.1E-3 ) 478 532 breakup = selfcoll * ( phi_br + 1.0 ) 479 533 ELSE … … 481 535 ENDIF 482 536 483 selfcoll = MAX( breakup - selfcoll, -nr(k,j,i) / ( dt_3d * & 484 weight_substep(intermediate_timestep_count) ) ) 485 ! 486 !-- Tendency for nr: 487 tend_nr(k,j,i) = tend_nr(k,j,i) + selfcoll 537 selfcoll = MAX( breakup - selfcoll, -nr_1d(k) / dt_micro ) 538 nr_1d(k) = nr_1d(k) + selfcoll * dt_micro 488 539 489 540 ENDIF 490 491 541 ENDDO 492 542 … … 502 552 USE cloud_parameters 503 553 USE constants 504 USE indices 505 USE control_parameters 506 USE statistics 507 508 IMPLICIT NONE 509 510 INTEGER :: i, j, k 511 REAL :: evap, alpha, e_s, q_s, t_l, sat, temp, g_evap, f_vent, & 512 mu_r_2, mu_r_5d2, nr_0 513 514 515 DO k = nzb_2d(j,i)+1, nzt 516 517 IF ( qr(k,j,i) > eps_sb ) THEN 554 USE control_parameters 555 USE indices 556 557 IMPLICIT NONE 558 559 INTEGER :: i, j, k 560 REAL :: alpha, dr, e_s, evap, evap_nr, f_vent, g_evap, lambda_r, & 561 mu_r, mu_r_2, mu_r_5d2, nr_0, q_s, sat, t_l, temp, xr 562 563 DO k = nzb_s_inner(j,i)+1, nzt 564 IF ( qr_1d(k) > eps_sb ) THEN 518 565 ! 519 566 !-- Actual liquid water temperature: 520 t_l = t_d_pt(k) * pt (k,j,i)567 t_l = t_d_pt(k) * pt_1d(k) 521 568 ! 522 569 !-- Saturation vapor pressure at t_l: … … 526 573 q_s = 0.622 * e_s / ( hyp(k) - 0.378 * e_s ) 527 574 alpha = 0.622 * l_d_r * l_d_cp / ( t_l * t_l ) 528 q_s = q_s * ( 1.0 + alpha * q (k,j,i) ) / ( 1.0 + alpha * q_s )575 q_s = q_s * ( 1.0 + alpha * q_1d(k) ) / ( 1.0 + alpha * q_s ) 529 576 ! 530 577 !-- Supersaturation: 531 sat = MIN( 0.0, ( q (k,j,i) - ql(k,j,i) ) / q_s - 1.0 )578 sat = MIN( 0.0, ( q_1d(k) - qr_1d(k) - qc_1d(k) ) / q_s - 1.0 ) 532 579 ! 533 580 !-- Actual temperature: 534 temp = t_l + l_d_cp * ql(k,j,i) 535 ! 536 !-- 537 g_evap = ( l_v / ( r_v * temp ) - 1.0 ) * l_v / & 538 ( thermal_conductivity_l * temp ) + rho_l * r_v * temp /& 539 ( diff_coeff_l * e_s ) 540 g_evap = 1.0 / g_evap 581 temp = t_l + l_d_cp * ( qc_1d(k) + qr_1d(k) ) 582 583 g_evap = 1.0 / ( ( l_v / ( r_v * temp ) - 1.0 ) * l_v / & 584 ( thermal_conductivity_l * temp ) + r_v * temp / & 585 ( diff_coeff_l * e_s ) ) 586 ! 587 !-- Mean weight of rain drops 588 xr = hyrho(k) * qr_1d(k) / nr_1d(k) 589 ! 590 !-- Weight averaged diameter of rain drops: 591 dr = ( xr * dpirho_l )**( 1.0 / 3.0 ) 541 592 ! 542 593 !-- Compute ventilation factor and intercept parameter 543 594 !-- (Seifert and Beheng, 2006; Seifert, 2008): 544 595 IF ( ventilation_effect ) THEN 545 mu_r_2 = mu_r(k) + 2.0 546 mu_r_5d2 = mu_r(k) + 2.5 596 ! 597 !-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005; 598 !-- Stevens and Seifert, 2008): 599 mu_r = 10.0 * ( 1.0 + TANH( 1.2E3 * ( dr - 1.4E-3 ) ) ) 600 ! 601 !-- Slope parameter of gamma distribution (Seifert, 2008): 602 lambda_r = ( ( mu_r + 3.0 ) * ( mu_r + 2.0 ) * & 603 ( mu_r + 1.0 ) )**( 1.0 / 3.0 ) / dr 604 605 mu_r_2 = mu_r + 2.0 606 mu_r_5d2 = mu_r + 2.5 547 607 f_vent = a_vent * gamm( mu_r_2 ) * & 548 lambda_r (k)**( -mu_r_2 ) +&608 lambda_r**( -mu_r_2 ) + & 549 609 b_vent * schmidt_p_1d3 * & 550 610 SQRT( a_term / kin_vis_air ) * gamm( mu_r_5d2 ) * & 551 lambda_r (k)**( -mu_r_5d2 ) *&611 lambda_r**( -mu_r_5d2 ) * & 552 612 ( 1.0 - 0.5 * ( b_term / a_term ) * & 553 ( lambda_r (k) /&554 ( c_term + lambda_r (k) ) )**mu_r_5d2 -&613 ( lambda_r / & 614 ( c_term + lambda_r ) )**mu_r_5d2 - & 555 615 0.125 * ( b_term / a_term )**2 * & 556 ( lambda_r (k) /&557 ( 2.0 * c_term + lambda_r (k) ) )**mu_r_5d2 -&616 ( lambda_r / & 617 ( 2.0 * c_term + lambda_r ) )**mu_r_5d2 - & 558 618 0.0625 * ( b_term / a_term )**3 * & 559 ( lambda_r (k) /&560 ( 3.0 * c_term + lambda_r (k) ) )**mu_r_5d2 -&619 ( lambda_r / & 620 ( 3.0 * c_term + lambda_r ) )**mu_r_5d2 - & 561 621 0.0390625 * ( b_term / a_term )**4 * & 562 ( lambda_r (k) /&563 ( 4.0 * c_term + lambda_r (k)) )**mu_r_5d2 )564 nr_0 = nr (k,j,i) * lambda_r(k)**( mu_r(k) + 1.0 ) /&565 gamm( mu_r (k)+ 1.0 )622 ( lambda_r / & 623 ( 4.0 * c_term + lambda_r ) )**mu_r_5d2 ) 624 nr_0 = nr_1d(k) * lambda_r**( mu_r + 1.0 ) / & 625 gamm( mu_r + 1.0 ) 566 626 ELSE 567 627 f_vent = 1.0 568 nr_0 = nr (k,j,i) * dr(k)628 nr_0 = nr_1d(k) * dr 569 629 ENDIF 570 630 ! … … 572 632 evap = 2.0 * pi * nr_0 * g_evap * f_vent * sat / & 573 633 hyrho(k) 574 evap = MAX( evap, -qr(k,j,i) / ( dt_3d * & 575 weight_substep(intermediate_timestep_count) ) ) 576 ! 577 !-- Tendencies for q, qr, nr, pt: 578 tend_qr(k,j,i) = tend_qr(k,j,i) + evap 579 tend_q(k,j,i) = tend_q(k,j,i) - evap 580 tend_nr(k,j,i) = tend_nr(k,j,i) + c_evap * evap / xr(k) * hyrho(k) 581 tend_pt(k,j,i) = tend_pt(k,j,i) + evap * l_d_cp * pt_d_t(k) 582 634 635 evap = MAX( evap, -qr_1d(k) / dt_micro ) 636 evap_nr = MAX( c_evap * evap / xr * hyrho(k), & 637 -nr_1d(k) / dt_micro ) 638 639 qr_1d(k) = qr_1d(k) + evap * dt_micro 640 nr_1d(k) = nr_1d(k) + evap_nr * dt_micro 583 641 ENDIF 584 642 … … 593 651 USE cloud_parameters 594 652 USE constants 595 USE indices596 USE control_parameters653 USE control_parameters 654 USE indices 597 655 598 656 IMPLICIT NONE 599 657 600 658 INTEGER :: i, j, k 601 REAL :: sed_q_const, sigma_gc = 1.3, k_st = 1.2E8 659 REAL :: sed_qc_const 660 661 REAL, DIMENSION(nzb:nzt+1) :: sed_qc 602 662 603 663 ! 604 664 !-- Sedimentation of cloud droplets (Heus et al., 2010): 605 sed_q _const = k_st * ( 3.0 / ( 4.0 * pi * rho_l ))**( 2.0 / 3.0 ) *&665 sed_qc_const = k_st * ( 3.0 / ( 4.0 * pi * rho_l ))**( 2.0 / 3.0 ) * & 606 666 EXP( 5.0 * LOG( sigma_gc )**2 ) 607 667 608 sed_q = 0.0 609 610 DO k = nzb_2d(j,i)+1, nzt 611 IF ( ql(k,j,i) > 0.0 ) THEN 612 sed_q(k) = sed_q_const * nc**( -2.0 / 3.0 ) * & 613 ( ql(k,j,i) * hyrho(k) )**( 5.0 / 3.0 ) 668 sed_qc(nzt+1) = 0.0 669 670 DO k = nzt, nzb_s_inner(j,i)+1, -1 671 IF ( qc_1d(k) > eps_sb ) THEN 672 sed_qc(k) = sed_qc_const * nc_1d(k)**( -2.0 / 3.0 ) * & 673 ( qc_1d(k) * hyrho(k) )**( 5.0 / 3.0 ) 674 ELSE 675 sed_qc(k) = 0.0 614 676 ENDIF 615 ENDDO 616 ! 617 !-- Tendency for q, pt: 618 DO k = nzb_2d(j,i)+1, nzt 619 tend_q(k,j,i) = tend_q(k,j,i) + ( sed_q(k+1) - sed_q(k) ) * & 620 ddzu(k+1) / hyrho(k) 621 tend_pt(k,j,i) = tend_pt(k,j,i) - ( sed_q(k+1) - sed_q(k) ) * & 622 ddzu(k+1) / hyrho(k) * l_d_cp * pt_d_t(k) 677 678 sed_qc(k) = MIN( sed_qc(k), hyrho(k) * dzu(k+1) * q_1d(k) / & 679 dt_micro + sed_qc(k+1) ) 680 681 q_1d(k) = q_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / & 682 hyrho(k) * dt_micro 683 qc_1d(k) = qc_1d(k) + ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / & 684 hyrho(k) * dt_micro 685 pt_1d(k) = pt_1d(k) - ( sed_qc(k+1) - sed_qc(k) ) * ddzu(k+1) / & 686 hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro 687 623 688 ENDDO 624 689 … … 631 696 USE cloud_parameters 632 697 USE constants 633 USE indices634 USE control_parameters698 USE control_parameters 699 USE indices 635 700 USE statistics 636 701 … … 638 703 639 704 INTEGER :: i, j, k, k_run 640 REAL :: c_run, d_max, d_mean, d_min, dt_sedi, flux, z_run 641 642 REAL, DIMENSION(nzb:nzt) :: c_nr, c_qr, nr_slope, qr_slope, w_nr, w_qr 643 705 REAL :: c_run, d_max, d_mean, d_min, dr, dt_sedi, flux, lambda_r, & 706 mu_r, z_run 707 708 REAL, DIMENSION(nzb:nzt+1) :: c_nr, c_qr, d_nr, d_qr, nr_slope, & 709 qr_slope, sed_nr, sed_qr, w_nr, w_qr 644 710 ! 645 711 !-- Computation of sedimentation flux. Implementation according to Stevens 646 712 !-- and Seifert (2008). 647 713 IF ( intermediate_timestep_count == 1 ) prr(:,j,i) = 0.0 648 649 dt_sedi = dt_3d * weight_substep(intermediate_timestep_count)650 651 w_nr = 0.0652 w_qr = 0.0653 714 ! 654 715 !-- Compute velocities 655 716 DO k = nzb_s_inner(j,i)+1, nzt 656 IF ( qr(k,j,i) > eps_sb ) THEN 717 IF ( qr_1d(k) > eps_sb ) THEN 718 ! 719 !-- Weight averaged diameter of rain drops: 720 dr = ( hyrho(k) * qr_1d(k) / nr_1d(k) * dpirho_l )**( 1.0 / 3.0 ) 721 ! 722 !-- Shape parameter of gamma distribution (Milbrandt and Yau, 2005; 723 !-- Stevens and Seifert, 2008): 724 mu_r = 10.0 * ( 1.0 + TANH( 1.2E3 * ( dr - 1.4E-3 ) ) ) 725 ! 726 !-- Slope parameter of gamma distribution (Seifert, 2008): 727 lambda_r = ( ( mu_r + 3.0 ) * ( mu_r + 2.0 ) * & 728 ( mu_r + 1.0 ) )**( 1.0 / 3.0 ) / dr 729 657 730 w_nr(k) = MAX( 0.1, MIN( 20.0, a_term - b_term * ( 1.0 + & 658 c_term / lambda_r (k) )**( -1.0 * ( mu_r(k)+ 1.0 ) ) ) )731 c_term / lambda_r )**( -1.0 * ( mu_r + 1.0 ) ) ) ) 659 732 w_qr(k) = MAX( 0.1, MIN( 20.0, a_term - b_term * ( 1.0 + & 660 c_term / lambda_r (k) )**( -1.0 * ( mu_r(k)+ 4.0 ) ) ) )733 c_term / lambda_r )**( -1.0 * ( mu_r + 4.0 ) ) ) ) 661 734 ELSE 662 735 w_nr(k) = 0.0 … … 666 739 ! 667 740 !-- Adjust boundary values 668 w_nr(nzb_ 2d(j,i)) = w_nr(nzb_2d(j,i)+1)669 w_qr(nzb_ 2d(j,i)) = w_qr(nzb_2d(j,i)+1)670 w_nr(nzt ) = w_nr(nzt-1)671 w_qr(nzt ) = w_qr(nzt-1)741 w_nr(nzb_s_inner(j,i)) = w_nr(nzb_s_inner(j,i)+1) 742 w_qr(nzb_s_inner(j,i)) = w_qr(nzb_s_inner(j,i)+1) 743 w_nr(nzt+1) = 0.0 744 w_qr(nzt+1) = 0.0 672 745 ! 673 746 !-- Compute Courant number 674 DO k = nzb_s_inner(j,i)+1, nzt -1747 DO k = nzb_s_inner(j,i)+1, nzt 675 748 c_nr(k) = 0.25 * ( w_nr(k-1) + 2.0 * w_nr(k) + w_nr(k+1) ) * & 676 dt_ sedi* ddzu(k)749 dt_micro * ddzu(k) 677 750 c_qr(k) = 0.25 * ( w_qr(k-1) + 2.0 * w_qr(k) + w_qr(k+1) ) * & 678 dt_ sedi* ddzu(k)679 ENDDO 751 dt_micro * ddzu(k) 752 ENDDO 680 753 ! 681 754 !-- Limit slopes with monotonized centered (MC) limiter (van Leer, 1977): 682 755 IF ( limiter_sedimentation ) THEN 683 756 684 qr(nzb_s_inner(j,i),j,i) = 0.0 685 nr(nzb_s_inner(j,i),j,i) = 0.0 686 qr(nzt,j,i) = 0.0 687 nr(nzt,j,i) = 0.0 688 689 DO k = nzb_s_inner(j,i)+1, nzt-1 690 d_mean = 0.5 * ( qr(k+1,j,i) + qr(k-1,j,i) ) 691 d_min = qr(k,j,i) - MIN( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) 692 d_max = MAX( qr(k+1,j,i), qr(k,j,i), qr(k-1,j,i) ) - qr(k,j,i) 757 DO k = nzb_s_inner(j,i)+1, nzt 758 d_mean = 0.5 * ( qr_1d(k+1) + qr_1d(k-1) ) 759 d_min = qr_1d(k) - MIN( qr_1d(k+1), qr_1d(k), qr_1d(k-1) ) 760 d_max = MAX( qr_1d(k+1), qr_1d(k), qr_1d(k-1) ) - qr_1d(k) 693 761 694 762 qr_slope(k) = SIGN(1.0, d_mean) * MIN ( 2.0 * d_min, 2.0 * d_max, & 695 763 ABS( d_mean ) ) 696 764 697 d_mean = 0.5 * ( nr (k+1,j,i) + nr(k-1,j,i) )698 d_min = nr (k,j,i) - MIN( nr(k+1,j,i), nr(k,j,i), nr(k-1,j,i) )699 d_max = MAX( nr (k+1,j,i), nr(k,j,i), nr(k-1,j,i) ) - nr(k,j,i)765 d_mean = 0.5 * ( nr_1d(k+1) + nr_1d(k-1) ) 766 d_min = nr_1d(k) - MIN( nr_1d(k+1), nr_1d(k), nr_1d(k-1) ) 767 d_max = MAX( nr_1d(k+1), nr_1d(k), nr_1d(k-1) ) - nr_1d(k) 700 768 701 769 nr_slope(k) = SIGN(1.0, d_mean) * MIN ( 2.0 * d_min, 2.0 * d_max, & … … 709 777 710 778 ENDIF 779 780 sed_nr(nzt+1) = 0.0 781 sed_qr(nzt+1) = 0.0 711 782 ! 712 783 !-- Compute sedimentation flux 713 DO k = nzt -2, nzb_s_inner(j,i)+1, -1784 DO k = nzt, nzb_s_inner(j,i)+1, -1 714 785 ! 715 786 !-- Sum up all rain drop number densities which contribute to the flux … … 719 790 k_run = k 720 791 c_run = MIN( 1.0, c_nr(k) ) 721 DO WHILE ( c_run > 0.0 .AND. k_run <= nzt -1)792 DO WHILE ( c_run > 0.0 .AND. k_run <= nzt ) 722 793 flux = flux + hyrho(k_run) * & 723 ( nr (k_run,j,i) + nr_slope(k_run) * ( 1.0 - c_run ) *&794 ( nr_1d(k_run) + nr_slope(k_run) * ( 1.0 - c_run ) * & 724 795 0.5 ) * c_run * dzu(k_run) 725 796 z_run = z_run + dzu(k_run) … … 731 802 !-- available 732 803 flux = MIN( flux, & 733 hyrho(k) * dzu(k ) * nr(k,j,i) + sed_nr(k+1) * dt_sedi)734 735 sed_nr(k) = flux / dt_sedi736 tend_nr(k,j,i) = tend_nr(k,j,i) + ( sed_nr(k+1) - sed_nr(k) ) *&737 ddzu(k+1) / hyrho(k)804 hyrho(k) * dzu(k+1) * nr_1d(k) + sed_nr(k+1) * dt_micro ) 805 806 sed_nr(k) = flux / dt_micro 807 nr_1d(k) = nr_1d(k) + ( sed_nr(k+1) - sed_nr(k) ) * ddzu(k+1) / & 808 hyrho(k) * dt_micro 738 809 ! 739 810 !-- Sum up all rain water content which contributes to the flux … … 747 818 748 819 flux = flux + hyrho(k_run) * & 749 ( qr (k_run,j,i) + qr_slope(k_run) * ( 1.0 - c_run ) * &820 ( qr_1d(k_run) + qr_slope(k_run) * ( 1.0 - c_run ) * & 750 821 0.5 ) * c_run * dzu(k_run) 751 822 z_run = z_run + dzu(k_run) … … 757 828 !-- It is not allowed to sediment more rain water content than available 758 829 flux = MIN( flux, & 759 hyrho(k) * dzu(k) * qr(k,j,i) + sed_qr(k+1) * dt_sedi ) 760 761 sed_qr(k) = flux / dt_sedi 762 tend_qr(k,j,i) = tend_qr(k,j,i) + ( sed_qr(k+1) - sed_qr(k) ) * & 763 ddzu(k+1) / hyrho(k) 830 hyrho(k) * dzu(k) * qr_1d(k) + sed_qr(k+1) * dt_micro ) 831 832 sed_qr(k) = flux / dt_micro 833 834 qr_1d(k) = qr_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / & 835 hyrho(k) * dt_micro 836 q_1d(k) = q_1d(k) + ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / & 837 hyrho(k) * dt_micro 838 pt_1d(k) = pt_1d(k) - ( sed_qr(k+1) - sed_qr(k) ) * ddzu(k+1) / & 839 hyrho(k) * l_d_cp * pt_d_t(k) * dt_micro 764 840 ! 765 841 !-- Compute the rain rate 766 842 prr(k,j,i) = prr(k,j,i) + sed_qr(k) / hyrho(k) * & 767 weight_substep(intermediate_timestep_count) 768 ENDDO 843 weight_substep(intermediate_timestep_count) 844 ENDDO 845 769 846 ! 770 847 !-- Precipitation amount … … 774 851 775 852 precipitation_amount(j,i) = precipitation_amount(j,i) + & 776 prr(nzb_ 2d(j,i)+1,j,i) * &777 hyrho(nzb_ 2d(j,i)+1) * dt_3d853 prr(nzb_s_inner(j,i)+1,j,i) * & 854 hyrho(nzb_s_inner(j,i)+1) * dt_3d 778 855 ENDIF 779 856
Note: See TracChangeset
for help on using the changeset viewer.
