Changeset 5406 for LMDZ6/trunk/libf/phylmd/lmdz_lscp_condensation.f90
- Timestamp:
- Dec 13, 2024, 11:40:07 AM (2 days ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/lmdz_lscp_condensation.f90
r5396 r5406 8 8 9 9 CONTAINS 10 11 10 12 11 !********************************************************************************** … … 122 121 USE lmdz_lscp_ini, ONLY: lunout 123 122 124 USE lmdz_lscp_ini, ONLY: depo_coef_cirrus, capa_cond_cirrus, &123 USE lmdz_lscp_ini, ONLY: depo_coef_cirrus, capa_cond_cirrus, std_subl_pdf_lscp, & 125 124 mu_subl_pdf_lscp, beta_pdf_lscp, temp_thresh_pdf_lscp, & 126 rhlmid_pdf_lscp, k0_pdf_lscp, kappa_pdf_lscp, rhl0_pdf_lscp, &125 std100_pdf_lscp, k0_pdf_lscp, kappa_pdf_lscp, & 127 126 coef_mixing_lscp, coef_shear_lscp, & 128 127 chi_mixing_lscp, rho_ice … … 151 150 REAL, INTENT(IN) , DIMENSION(klon) :: qsat ! saturation specific humidity [kg/kg] 152 151 REAL, INTENT(IN) , DIMENSION(klon) :: gamma_cond ! condensation threshold w.r.t. qsat [-] 153 REAL, INTENT(IN OUT), DIMENSION(klon) :: ratqs ! ratio between the variance of the total water distribution and its average [-]152 REAL, INTENT(IN) , DIMENSION(klon) :: ratqs ! ratio between the variance of the total water distribution and its average [-] 154 153 LOGICAL, INTENT(IN) , DIMENSION(klon) :: keepgoing ! .TRUE. if a new condensation loop should be computed 155 154 ! … … 211 210 ! for unadjusted clouds 212 211 REAL :: qvapincld, qvapincld_new 213 REAL :: qiceincld, qice_ratio 214 REAL :: pres_sat, kappa 215 REAL :: air_thermal_conduct, water_vapor_diff 216 REAL :: iwc 217 REAL :: iwc_log_inf100, iwc_log_sup100 218 REAL :: iwc_inf100, alpha_inf100, coef_inf100 219 REAL :: mu_sup100, sigma_sup100, coef_sup100 220 REAL :: Dm_ice, rm_ice 212 REAL :: qiceincld 221 213 ! 222 214 ! for sublimation … … 309 301 310 302 pdf_std = ratqs(i) * qtot(i) 311 pdf_k = -SQRT( LOG( 1. + (pdf_std / qtot(i))**2 .) )303 pdf_k = -SQRT( LOG( 1. + (pdf_std / qtot(i))**2 ) ) 312 304 pdf_delta = LOG( qtot(i) / ( gamma_cond(i) * qsat(i) ) ) 313 305 pdf_a = pdf_delta / ( pdf_k * SQRT(2.) ) … … 342 334 !--and the condensation process is slightly adapted 343 335 !--This can happen only if ok_weibull_warm_clouds = .TRUE. 344 cldfra(i) = 0. 345 qvc(i) = 0. 346 qcld(i) = 0. 336 ! AB WARM CLOUD 337 !cldfra(i) = 0. 338 !qcld(i) = 0. 339 !qvc(i) = 0. 340 cldfra(i) = MAX(0., MIN(1., cf_seri(i))) 341 qcld(i) = MAX(0., MIN(qtot(i), ql_seri(i) + qi_seri(i) + rvc_seri(i) * qtot(i))) 342 qvc(i) = MAX(0., MIN(qcld(i), rvc_seri(i) * qtot(i))) 347 343 ok_warm_cloud = .TRUE. 348 344 ELSE … … 384 380 !--Cell dry air mass [kg] 385 381 !M_cell = rhodz * cell_area(i) 386 387 388 IF ( ok_unadjusted_clouds ) THEN389 !--If ok_unadjusted_clouds is set to TRUE, then the saturation adjustment390 !--hypothesis is lost, and the vapor in the cloud is purely prognostic.391 !392 !--The deposition equation is393 !-- dmi/dt = alpha*4pi*C*Svi / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) )394 !--from Lohmann et al. (2016), where395 !--alpha is the deposition coefficient [-]396 !--mi is the mass of one ice crystal [kg]397 !--C is the capacitance of an ice crystal [m]398 !--Svi is the supersaturation ratio equal to (qvc - qsat)/qsat [-]399 !--R_v is the specific gas constant for humid air [J/kg/K]400 !--T is the temperature [K]401 !--esi is the saturation pressure w.r.t. ice [Pa]402 !--Dv is the diffusivity of water vapor [m2/s]403 !--Ls is the specific latent heat of sublimation [J/kg/K]404 !--ka is the thermal conductivity of dry air [J/m/s/K]405 !406 !--alpha is a coefficient to take into account the fact that during deposition, a water407 !--molecule cannot join the crystal from everywhere, it must do so that the crystal stays408 !--coherent (with the same structure). It has no impact for sublimation.409 !--We fix alpha = depo_coef_cirrus (=0.5 by default following Lohmann et al. (2016))410 !--during deposition, and alpha = 1. during sublimation.411 !--The capacitance of the ice crystals is proportional to a parameter capa_cond_cirrus412 !-- C = capa_cond_cirrus * rm_ice413 !414 !--We have qice = Nice * mi, where Nice is the ice crystal415 !--number concentration per kg of moist air416 !--HYPOTHESIS 1: the ice crystals are spherical, therefore417 !-- mi = 4/3 * pi * rm_ice**3 * rho_ice418 !--HYPOTHESIS 2: the ice crystals are monodisperse with the419 !--initial radius rm_ice_0.420 !--NB. this is notably different than the assumption421 !--of a distributed qice in the cloud made in the sublimation process;422 !--should it be consistent?423 !424 !--As the deposition process does not create new ice crystals,425 !--and because we assume a same rm_ice value for all crystals426 !--therefore the sublimation process does not destroy ice crystals427 !--(or, in a limit case, it destroys all ice crystals), then428 !--Nice is a constant during the sublimation/deposition process.429 !-- dmi = dqi, et Nice = qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )430 !431 !--The deposition equation then reads:432 !-- dqi/dt = alpha*4pi*capa_cond_cirrus*rm_ice*(qvc-qsat)/qsat / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) ) * Nice433 !-- dqi/dt = alpha*4pi*capa_cond_cirrus* (qi / qi_0)**(1/3) *rm_ice_0*(qvc-qsat)/qsat &434 !-- / ( R_v*T/esi/Dv + Ls/ka/T * (Ls*R_v/T - 1) ) &435 !-- * qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice )436 !-- dqi/dt = qi**(1/3) * (qvc - qsat) * qi_0**(2/3) &437 !-- *alpha/qsat*capa_cond_cirrus/ (R_v*T/esi/Dv + Ls/ka/T*(Ls*R_v/T - 1)) / ( 1/3 rm_ice_0**2 rho_ice )438 !--and we have439 !-- dqvc/dt = - qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2440 !-- dqi/dt = qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2441 !--where kappa = 1/3*rho_ice/capa_cond_cirrus*qsat*(R_v*T/esi/Dv + Ls/ka/T*(Ls/R_v/T - 1))442 !443 !--This system of equations can be resolved with an exact444 !--explicit numerical integration, having one variable resolved445 !--explicitly, the other exactly. The exactly resolved variable is446 !--the one decreasing, so it is qvc if the process is447 !--condensation, qi if it is sublimation.448 !449 !--kappa is computed as an initialisation constant, as it depends only450 !--on temperature and other pre-computed values451 pres_sat = qsat(i) / ( EPS_W + ( 1. - EPS_W ) * qsat(i) ) * pplay(i)452 !--This formula for air thermal conductivity comes from Beard and Pruppacher (1971)453 air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184454 !--This formula for water vapor diffusivity comes from Hall and Pruppacher (1976)455 water_vapor_diff = 0.211 * ( temp(i) / RTT )**1.94 * ( 101325. / pplay(i) ) * 1.e-4456 kappa = 1. / 3. * rho_ice / capa_cond_cirrus * qsat(i) &457 * ( RV * temp(i) / water_vapor_diff / pres_sat &458 + RLSTT / air_thermal_conduct / temp(i) * ( RLSTT / RV / temp(i) - 1. ) )459 !--NB. the greater kappa, the lower the efficiency of the deposition/sublimation process460 ENDIF461 382 462 383 … … 491 412 qiceincld = ( qcld(i) / cldfra(i) - qvapincld ) 492 413 493 IF ( ok_unadjusted_clouds ) THEN 494 !--Here, the initial vapor in the cloud is qvapincld, and we compute 495 !--the new vapor qvapincld_new 496 497 !--rm_ice formula from McFarquhar and Heymsfield (1997) 498 iwc = qiceincld * rho * 1e3 499 iwc_inf100 = MIN(iwc, 0.252 * iwc**0.837) 500 iwc_log_inf100 = LOG10( MAX(eps, iwc_inf100) ) 501 iwc_log_sup100 = LOG10( MAX(eps, iwc - iwc_inf100) ) 502 503 alpha_inf100 = - 4.99E-3 - 0.0494 * iwc_log_inf100 504 coef_inf100 = iwc_inf100 * alpha_inf100**3. / 120. 505 506 mu_sup100 = ( 5.2 + 0.0013 * ( temp(i) - RTT ) ) & 507 + ( 0.026 - 1.2E-3 * ( temp(i) - RTT ) ) * iwc_log_sup100 508 sigma_sup100 = ( 0.47 + 2.1E-3 * ( temp(i) - RTT ) ) & 509 + ( 0.018 - 2.1E-4 * ( temp(i) - RTT ) ) * iwc_log_sup100 510 coef_sup100 = ( iwc - iwc_inf100 ) / EXP( 3 * mu_sup100 + 4.5 * sigma_sup100**2. ) 511 512 Dm_ice = ( 2. / alpha_inf100 * coef_inf100 + EXP( mu_sup100 + 0.5 * sigma_sup100**2. ) & 513 * coef_sup100 ) / ( coef_inf100 + coef_sup100 ) 514 rm_ice = Dm_ice / 2. * 1.E-6 515 516 IF ( qvapincld .GE. qsat(i) ) THEN 517 !--If the cloud is initially supersaturated 518 !--Exact explicit integration (qvc exact, qice explicit) 519 qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) & 520 * EXP( - depo_coef_cirrus * dtime * qiceincld / kappa / rm_ice**2. ) 521 ELSE 522 !--If the cloud is initially subsaturated 523 !--Exact explicit integration (qice exact, qvc explicit) 524 !--The barrier is set so that the resulting vapor in cloud 525 !--cannot be greater than qsat 526 !--qice_ratio is the ratio between the new ice content and 527 !--the old one, it is comprised between 0 and 1 528 qice_ratio = ( 1. - 2. / 3. / kappa / rm_ice**2. * dtime * ( qsat(i) - qvapincld ) ) 529 530 IF ( qice_ratio .LT. 0. ) THEN 531 !--If all the ice has been sublimated, we sublimate 532 !--completely the cloud and do not activate the sublimation 533 !--process 534 !--Tendencies and diagnostics 535 dcf_sub(i) = - cldfra(i) 536 dqvc_sub(i) = - qvc(i) 537 dqi_sub(i) = - ( qcld(i) - qvc(i) ) 538 539 cldfra(i) = 0. 540 qcld(i) = 0. 541 qvc(i) = 0. 542 543 !--The new vapor in cloud is set to 0 so that the 544 !--sublimation process does not activate 545 qvapincld_new = 0. 546 ELSE 547 !--Else, the sublimation process is activated with the 548 !--diagnosed new cloud water vapor 549 !--The new vapor in the cloud is increased with the 550 !--sublimated ice 551 qvapincld_new = qvapincld + qiceincld * ( 1. - qice_ratio**(3./2.) ) 552 !--The new vapor in the cloud cannot be greater than qsat 553 qvapincld_new = MIN(qvapincld_new, qsat(i)) 554 ENDIF ! qice_ratio .LT. 0. 555 ENDIF ! qvapincld .GT. qsat(i) 414 ! AB WARM CLOUD 415 !IF ( ok_unadjusted_clouds ) THEN 416 IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN 417 CALL deposition_sublimation(qvapincld, qiceincld, temp(i), qsat(i), & 418 pplay(i), dtime, qvapincld_new) 419 IF ( qvapincld_new .EQ. 0. ) THEN 420 !--If all the ice has been sublimated, we sublimate 421 !--completely the cloud and do not activate the sublimation 422 !--process 423 !--Tendencies and diagnostics 424 dcf_sub(i) = - cldfra(i) 425 dqvc_sub(i) = - qvc(i) 426 dqi_sub(i) = - ( qcld(i) - qvc(i) ) 427 428 cldfra(i) = 0. 429 qcld(i) = 0. 430 qvc(i) = 0. 431 ENDIF 556 432 ELSE 557 433 !--We keep the saturation adjustment hypothesis, and the vapor in the … … 576 452 577 453 !--If the vapor in cloud is below vapor needed for the cloud to survive 578 IF ( qvapincld .LT. qvapincld_new) THEN454 IF ( ( qvapincld .LT. qvapincld_new ) .OR. ( iflag_cloud_sublim_pdf .GE. 4 ) ) THEN 579 455 !--Sublimation of the subsaturated cloud 580 456 !--iflag_cloud_sublim_pdf selects the PDF of the ice water content … … 589 465 590 466 dcf_sub(i) = - cldfra(i) * ( qvapincld_new - qvapincld ) * pdf_e1 591 dqt_sub = - cldfra(i) * ( qvapincld_new**2 . - qvapincld**2.) &467 dqt_sub = - cldfra(i) * ( qvapincld_new**2 - qvapincld**2 ) & 592 468 * pdf_e1 / 2. 593 469 … … 601 477 + qvapincld - qvapincld_new * pdf_e1 ) 602 478 603 ELSEIF ( iflag_cloud_sublim_pdf . GE. 3 ) THEN479 ELSEIF ( iflag_cloud_sublim_pdf .EQ. 3 ) THEN 604 480 !--Gamma distribution starting at qvapincld 605 481 pdf_alpha = ( mu_subl_pdf_lscp + 1. ) / qiceincld … … 610 486 dcf_sub(i) = - cldfra(i) * pdf_e1 611 487 dqt_sub = - cldfra(i) * ( pdf_e2 / pdf_alpha + qvapincld * pdf_e1 ) 488 489 ELSEIF ( iflag_cloud_sublim_pdf .EQ. 4 ) THEN 490 !--Normal distribution around qvapincld + qiceincld with width 491 !--constant in the RHi space 492 pdf_y = ( qvapincld_new - qvapincld - qiceincld ) & 493 / ( std_subl_pdf_lscp / 100. * qsat(i)) 494 pdf_e1 = 0.5 * ( 1. + ERF( pdf_y / SQRT(2.) ) ) 495 pdf_e2 = EXP( - pdf_y**2 / 2. ) / SQRT( 2. * RPI ) 496 497 dcf_sub(i) = - cldfra(i) * pdf_e1 498 dqt_sub = - cldfra(i) * ( ( qvapincld + qiceincld ) * pdf_e1 & 499 - std_subl_pdf_lscp / 100. * qsat(i) * pdf_e2 ) 500 612 501 ENDIF 613 502 … … 640 529 !--New PDF 641 530 rhl_clr = qclr / ( 1. - cldfra(i) ) / qsatl(i) * 100. 642 rhl_clr = MIN(rhl_clr, 2. * rhlmid_pdf_lscp)643 531 644 532 !--Calculation of the properties of the PDF … … 648 536 !--Coefficient for standard deviation: 649 537 !-- tuning coef * (clear sky area**0.25) * (function of temperature) 650 pdf_e1 = beta_pdf_lscp & 651 * ( ( 1. - cldfra(i) ) * cell_area(i) )**( 1. / 4. ) & 652 * MAX( temp(i) - temp_thresh_pdf_lscp, 0. ) 653 IF ( rhl_clr .GT. rhlmid_pdf_lscp ) THEN 654 pdf_std = pdf_e1 * ( 2. * rhlmid_pdf_lscp - rhl_clr ) / rhlmid_pdf_lscp 538 !pdf_e1 = beta_pdf_lscp & 539 ! * ( ( 1. - cldfra(i) ) * cell_area(i) )**( 1. / 4. ) & 540 ! * MAX( temp(i) - temp_thresh_pdf_lscp, 0. ) 541 !-- tuning coef * (cell area**0.25) * (function of temperature) 542 pdf_e1 = beta_pdf_lscp * ( ( 1. - cldfra(i) ) * cell_area(i) )**0.25 & 543 * MAX( temp(i) - temp_thresh_pdf_lscp, 0. ) 544 IF ( rhl_clr .GT. 50. ) THEN 545 pdf_std = ( pdf_e1 - std100_pdf_lscp ) * ( 100. - rhl_clr ) / 50. + std100_pdf_lscp 655 546 ELSE 656 pdf_std = pdf_e1 * rhl_clr / rhlmid_pdf_lscp547 pdf_std = pdf_e1 * rhl_clr / 50. 657 548 ENDIF 658 549 pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. ) 659 pdf_alpha = EXP( rhl_clr / rhl0_pdf_lscp ) * pdf_e3 550 pdf_alpha = EXP( rhl_clr / 100. ) * pdf_e3 551 pdf_alpha = MIN(10., pdf_alpha) 552 553 IF ( ok_warm_cloud ) THEN 554 !--If the statistical scheme is activated, the standard deviation is adapted 555 !--to depend on the pressure level. It is multiplied by ratqs, so that near the 556 !--surface std is almost 0, and upper than about 450 hPa the std is left untouched 557 pdf_std = pdf_std * ratqs(i) 558 ENDIF 660 559 661 560 pdf_e2 = GAMMA(1. + 1. / pdf_alpha) 662 pdf_scale = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2 .))561 pdf_scale = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2 )) 663 562 pdf_loc = rhl_clr - pdf_scale * pdf_e2 664 665 !--Diagnostics of ratqs666 ratqs(i) = pdf_std / ( qclr / ( 1. - cldfra(i) ) / qsatl(i) * 100. )667 563 668 564 !--Calculation of the newly condensed water and fraction (pronostic) … … 677 573 qt_cond = ( pdf_e3 + pdf_loc * cf_cond ) * qsatl(i) / 100. 678 574 679 IF ( ok_warm_cloud ) THEN 680 !--If the statistical scheme is activated, the calculated increase is equal 681 !--to the cloud fraction, we assume saturation adjustment, and the 682 !--condensation process stops 683 cldfra(i) = cf_cond 684 qcld(i) = qt_cond 685 qvc(i) = cldfra(i) * qsat(i) 686 687 ELSEIF ( cf_cond .GT. eps ) THEN 575 ! AB WARM CLOUD 576 !IF ( ok_warm_cloud ) THEN 577 ! !--If the statistical scheme is activated, the calculated increase is equal 578 ! !--to the cloud fraction, we assume saturation adjustment, and the 579 ! !--condensation process stops 580 ! cldfra(i) = cf_cond 581 ! qcld(i) = qt_cond 582 ! qvc(i) = cldfra(i) * qsat(i) 583 584 !ELSEIF ( cf_cond .GT. eps ) THEN 585 IF ( cf_cond .GT. eps ) THEN 688 586 689 587 dcf_con(i) = ( 1. - cldfra(i) ) * cf_cond … … 695 593 696 594 697 IF ( ok_unadjusted_clouds ) THEN 595 ! AB WARM CLOUD 596 !IF ( ok_unadjusted_clouds ) THEN 597 IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN 698 598 !--Here, the initial vapor in the cloud is gamma_cond*qsat, and we compute 699 599 !--the new vapor qvapincld. The timestep is divided by two because we do not … … 701 601 qvapincld = gamma_cond(i) * qsat(i) 702 602 qiceincld = dqt_con / dcf_con(i) - gamma_cond(i) * qsat(i) 703 704 !--rm_ice formula from McFarquhar and Heymsfield (1997) 705 iwc = qiceincld * rho * 1e3 706 iwc_inf100 = MIN(iwc, 0.252 * iwc**0.837) 707 iwc_log_inf100 = LOG10( MAX(eps, iwc_inf100) ) 708 iwc_log_sup100 = LOG10( MAX(eps, iwc - iwc_inf100) ) 709 710 alpha_inf100 = - 4.99E-3 - 0.0494 * iwc_log_inf100 711 coef_inf100 = iwc_inf100 * alpha_inf100**3. / 120. 712 713 mu_sup100 = ( 5.2 + 0.0013 * ( temp(i) - RTT ) ) & 714 + ( 0.026 - 1.2E-3 * ( temp(i) - RTT ) ) * iwc_log_sup100 715 sigma_sup100 = ( 0.47 + 2.1E-3 * ( temp(i) - RTT ) ) & 716 + ( 0.018 - 2.1E-4 * ( temp(i) - RTT ) ) * iwc_log_sup100 717 coef_sup100 = ( iwc - iwc_inf100 ) / EXP( 3. * mu_sup100 + 4.5 * sigma_sup100**2. ) 718 719 Dm_ice = ( 2. / alpha_inf100 * coef_inf100 + EXP( mu_sup100 + 0.5 * sigma_sup100**2. ) & 720 * coef_sup100 ) / ( coef_inf100 + coef_sup100 ) 721 rm_ice = Dm_ice / 2. * 1.E-6 722 !--As qvapincld is necessarily greater than qsat, we only 723 !--use the exact explicit formulation 724 !--Exact explicit version 725 qvapincld = qsat(i) + ( qvapincld - qsat(i) ) & 726 * EXP( - depo_coef_cirrus * dtime / 2. * qiceincld / kappa / rm_ice**2. ) 603 CALL deposition_sublimation(qvapincld, qiceincld, temp(i), qsat(i), & 604 pplay(i), dtime / 2., qvapincld_new) 605 qvapincld = qvapincld_new 727 606 ELSE 728 607 !--We keep the saturation adjustment hypothesis, and the vapor in the … … 735 614 dqi_con(i) = dqt_con - dqvc_con(i) 736 615 737 !-- Note that the tendencies are NOT added because they are738 !--added after the mixing process. In the following, the gridbox fraction is739 !-- 1. - dcf_con(i), and the total water in the gridbox is740 !-- qtot(i) - dqi_con(i) - dqvc_con(i)616 !--Add tendencies 617 cldfra(i) = MIN(1., cldfra(i) + dcf_con(i)) 618 qcld(i) = MIN(qtot(i), qcld(i) + dqt_con) 619 qvc(i) = MIN(qcld(i), qvc(i) + dqvc_con(i)) 741 620 742 621 ENDIF ! ok_warm_cloud, cf_cond .GT. eps 743 ENDIF ! ( 1. - cldfra(i) ) .GT. eps 744 745 !--If there is still clear sky, we diagnose the ISSR 746 !--We recalculte the PDF properties (after the condensation process) 747 IF ( ( ( 1. - dcf_con(i) - cldfra(i) ) .GT. eps ) .AND. .NOT. ok_warm_cloud ) THEN 748 !--Water quantity in the clear-sky + potential liquid cloud (gridbox average) 749 qclr = qtot(i) - dqi_con(i) - dqvc_con(i) - qcld(i) 750 751 !--New PDF 752 rhl_clr = qclr / ( 1. - dcf_con(i) - cldfra(i) ) / qsatl(i) * 100. 753 rhl_clr = MIN(rhl_clr, 2. * rhlmid_pdf_lscp) 754 755 !--Calculation of the properties of the PDF 756 !--Parameterization from IAGOS observations 757 !--pdf_e1 and pdf_e2 will be reused below 758 759 !--Coefficient for standard deviation: 760 !-- tuning coef * (clear sky area**0.25) * (function of temperature) 761 pdf_e1 = beta_pdf_lscp & 762 * ( ( 1. - dcf_con(i) - cldfra(i) ) * cell_area(i) )**( 1. / 4. ) & 763 * MAX( temp(i) - temp_thresh_pdf_lscp, 0. ) 764 IF ( rhl_clr .GT. rhlmid_pdf_lscp ) THEN 765 pdf_std = pdf_e1 * ( 2. * rhlmid_pdf_lscp - rhl_clr ) / rhlmid_pdf_lscp 766 ELSE 767 pdf_std = pdf_e1 * rhl_clr / rhlmid_pdf_lscp 768 ENDIF 769 pdf_e3 = k0_pdf_lscp + kappa_pdf_lscp * MAX( temp_nowater - temp(i), 0. ) 770 pdf_alpha = EXP( rhl_clr / rhl0_pdf_lscp ) * pdf_e3 771 772 pdf_e2 = GAMMA(1. + 1. / pdf_alpha) 773 pdf_scale = MAX(eps, pdf_std / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2. )) 774 pdf_loc = rhl_clr - pdf_scale * pdf_e2 775 622 776 623 !--We then calculate the part that is greater than qsat 777 !--and consider it supersaturated 778 624 !--and lower than gamma_cond * qsat, which is the ice supersaturated region 779 625 pdf_x = qsat(i) / qsatl(i) * 100. 780 626 pdf_y = ( MAX( pdf_x - pdf_loc, 0. ) / pdf_scale ) ** pdf_alpha 781 627 pdf_e3 = GAMMAINC ( 1. + 1. / pdf_alpha , pdf_y ) 782 628 pdf_e3 = pdf_scale * ( 1. - pdf_e3 ) * pdf_e2 783 issrfra(i) = EXP( - pdf_y ) * ( 1. - dcf_con(i) - cldfra(i) ) 784 qissr(i) = ( pdf_e3 * ( 1. - dcf_con(i) - cldfra(i) ) + pdf_loc * issrfra(i) ) * qsatl(i) / 100. 785 ENDIF 629 issrfra(i) = EXP( - pdf_y ) * ( 1. - cldfra(i) ) 630 qissr(i) = ( pdf_e3 * ( 1. - cldfra(i) ) + pdf_loc * issrfra(i) ) * qsatl(i) / 100. 631 632 issrfra(i) = issrfra(i) - dcf_con(i) 633 qissr(i) = qissr(i) - dqvc_con(i) - dqi_con(i) 634 635 ENDIF ! ( 1. - cldfra(i) ) .GT. eps 786 636 787 637 !--Calculation of the subsaturated clear sky fraction and water 788 subfra(i) = 1. - dcf_con(i) -cldfra(i) - issrfra(i)789 qsub(i) = qtot(i) - dqi_con(i) - dqvc_con(i) -qcld(i) - qissr(i)638 subfra(i) = 1. - cldfra(i) - issrfra(i) 639 qsub(i) = qtot(i) - qcld(i) - qissr(i) 790 640 791 641 … … 797 647 !--but does not cover the entire mesh. 798 648 ! 799 IF ( ( cldfra(i) .LT. ( 1. - dcf_con(i) - eps ) ) .AND. ( cldfra(i) .GT. eps ) & 800 .AND. .NOT. ok_warm_cloud ) THEN 649 ! AB WARM CLOUD 650 !IF ( ( cldfra(i) .LT. ( 1. - dcf_con(i) - eps ) ) .AND. ( cldfra(i) .GT. eps ) & 651 ! .AND. .NOT. ok_warm_cloud ) THEN 652 IF ( ( cldfra(i) .LT. ( 1. - eps ) ) .AND. ( cldfra(i) .GT. eps ) ) THEN 801 653 802 654 !--Initialisation … … 845 697 !-- * ( 3. * ( 1. + b / a ) - SQRT( (3. + b / a) * (1. + 3. * b / a) ) ) 846 698 !--and finally 847 N_cld_mix = coef_mixing_lscp * cldfra(i) * ( 1. - dcf_con(i) - cldfra(i) )**2.&848 * cell_area(i) * ( 1. - dcf_con(i) ) *bovera / RPI &849 / ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) )**2 .850 !--where coef_mix_lscp = ( alpha * norm_constant )**2 .699 N_cld_mix = coef_mixing_lscp * cldfra(i) * ( 1. - cldfra(i) )**2 & 700 * cell_area(i) * bovera / RPI & 701 / ( 3. * (1. + bovera) - SQRT( (3. + bovera) * (1. + 3. * bovera) ) )**2 702 !--where coef_mix_lscp = ( alpha * norm_constant )**2 851 703 !--N_cld_mix is the number of clouds in contact with clear sky, and can be non-integer 852 704 !--In particular, it is 0 if cldfra = 1 853 a_mix = SQRT( cell_area(i) * ( 1. - dcf_con(i) ) *cldfra(i) / bovera / N_cld_mix / RPI )705 a_mix = SQRT( cell_area(i) * cldfra(i) / bovera / N_cld_mix / RPI ) 854 706 855 707 !--The time required for turbulent diffusion to homogenize a region of size 856 !--L_mix is defined as (L_mix**2 ./tke_dissip)**(1./3.) (Pope, 2000; Field et al., 2014)708 !--L_mix is defined as (L_mix**2/tke_dissip)**(1./3.) (Pope, 2000; Field et al., 2014) 857 709 !--We compute L_mix and assume that the cloud is mixed over this length 858 L_mix = SQRT( dtime**3 .* pbl_eps(i) )710 L_mix = SQRT( dtime**3 * pbl_eps(i) ) 859 711 !--The mixing length cannot be greater than the semi-minor axis. In this case, 860 712 !--the entire cloud is mixed. … … 863 715 !--The fraction of clear sky mixed is 864 716 !-- N_cld_mix * ( (a + L_mix) * (b + L_mix) - a * b ) * RPI / cell_area 865 envfra_mix = N_cld_mix * RPI / cell_area(i) / ( 1. - dcf_con(i) )&866 * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2 .)717 envfra_mix = N_cld_mix * RPI / cell_area(i) & 718 * ( a_mix * ( 1. + bovera ) * L_mix + L_mix**2 ) 867 719 !--The fraction of cloudy sky mixed is 868 720 !-- N_cld_mix * ( a * b - (a - L_mix) * (b - L_mix) ) * RPI / cell_area 869 cldfra_mix = N_cld_mix * RPI / cell_area(i) / ( 1. - dcf_con(i) )&870 * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2 .)721 cldfra_mix = N_cld_mix * RPI / cell_area(i) & 722 * ( a_mix * ( 1. + bovera ) * L_mix - L_mix**2 ) 871 723 872 724 … … 878 730 !--The increase in size is 879 731 L_shear = coef_shear_lscp * shear(i) * dz * dtime 880 !--therefore, the total increase in fractionis732 !--therefore, the fraction of clear sky mixed is 881 733 !-- N_cld_mix * ( (a + L_shear) * b - a * b ) * RPI / 2. / cell_area 882 shear_fra = RPI * L_shear * a_mix * bovera / 2. * N_cld_mix & 883 / cell_area(i) / ( 1. - dcf_con(i) ) 734 !--and the fraction of cloud mixed is 735 !-- N_cld_mix * ( (a * b) - (a - L_shear) * b ) * RPI / 2. / cell_area 736 !--(note that they are equal) 737 shear_fra = RPI * L_shear * a_mix * bovera / 2. * N_cld_mix / cell_area(i) 884 738 !--and the environment and cloud mixed fractions are the same, 885 739 !--which we add to the previous calculated mixed fractions. … … 906 760 !--to this fraction (sigma_mix * cldfra_mix). 907 761 IF ( subfra(i) .GT. eps ) THEN 908 909 IF ( ok_unadjusted_clouds ) THEN 910 !--The subsaturated air is simply added to the cloud, 911 !--with the corresponding cloud fraction 912 !--If the cloud is too subsaturated, the sublimation process 913 !--activated in the following timestep will reduce the cloud 914 !--fraction 915 dcf_mix_sub = subfra_mix 916 dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i) 917 dqvc_mix_sub = dqt_mix_sub 918 762 !--We compute the total humidity in the mixed air, which 763 !--can be either sub- or supersaturated. 764 qvapinmix = ( qsub(i) * subfra_mix / subfra(i) & 765 + qcld(i) * cldfra_mix * sigma_mix / cldfra(i) ) & 766 / ( subfra_mix + cldfra_mix * sigma_mix ) 767 768 ! AB WARM CLOUD 769 !IF ( ok_unadjusted_clouds ) THEN 770 IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN 771 qiceincld = ( qcld(i) - qvc(i) ) * cldfra_mix * sigma_mix / cldfra(i) & 772 / ( subfra_mix + cldfra_mix * sigma_mix ) 773 CALL deposition_sublimation(qvapinmix, qiceincld, temp(i), qsat(i), & 774 pplay(i), dtime, qvapincld_new) 775 IF ( qvapincld_new .EQ. 0. ) THEN 776 !--If all the ice has been sublimated, we sublimate 777 !--completely the cloud and do not activate the sublimation 778 !--process 779 !--Tendencies and diagnostics 780 dcf_mix_sub = - sigma_mix * cldfra_mix 781 dqt_mix_sub = dcf_mix_sub * qcld(i) / cldfra(i) 782 dqvc_mix_sub = dcf_mix_sub * qvc(i) / cldfra(i) 783 ELSE 784 dcf_mix_sub = subfra_mix 785 dqt_mix_sub = dcf_mix_sub * qsub(i) / subfra(i) 786 dqvc_mix_sub = dcf_mix_sub * qvapincld_new 787 ENDIF 919 788 ELSE 920 !--We compute the total humidity in the mixed air, which921 !--can be either sub- or supersaturated.922 qvapinmix = ( qsub(i) * subfra_mix / subfra(i) &923 + qcld(i) * cldfra_mix * sigma_mix / cldfra(i) ) &924 / ( subfra_mix + cldfra_mix * sigma_mix )925 926 789 IF ( qvapinmix .GT. qsat(i) ) THEN 927 790 !--If the mixed air is supersaturated, we condense the subsaturated … … 944 807 IF ( issrfra(i) .GT. eps ) THEN 945 808 946 IF ( ok_unadjusted_clouds ) THEN 947 !--The ice supersaturated air is simply added to the 948 !--cloud, and supersaturated vapor will be deposited on the 949 !--cloud ice crystals by the deposition process in the 950 !--following timestep 809 ! AB WARM CLOUD 810 !IF ( ok_unadjusted_clouds ) THEN 811 IF ( ok_unadjusted_clouds .AND. .NOT. ok_warm_cloud ) THEN 812 qvapinmix = ( qissr(i) * issrfra_mix / issrfra(i) & 813 + qcld(i) * cldfra_mix * (1. - sigma_mix) / cldfra(i) ) & 814 / ( issrfra_mix + cldfra_mix * (1. - sigma_mix) ) 815 qiceincld = ( qcld(i) - qvc(i) ) * cldfra_mix * (1. - sigma_mix) / cldfra(i) & 816 / ( issrfra_mix + cldfra_mix * (1. - sigma_mix) ) 817 CALL deposition_sublimation(qvapinmix, qiceincld, temp(i), qsat(i), & 818 pplay(i), dtime, qvapincld_new) 951 819 dcf_mix_issr = issrfra_mix 952 820 dqt_mix_issr = dcf_mix_issr * qissr(i) / issrfra(i) 953 dqvc_mix_issr = d qt_mix_issr821 dqvc_mix_issr = dcf_mix_issr * qvapincld_new 954 822 ELSE 955 823 !--In this case, the additionnal vapor condenses … … 971 839 issrfra(i) = MAX(0., issrfra(i) - dcf_mix_issr) 972 840 qissr(i) = MAX(0., qissr(i) - dqt_mix_issr) 973 cldfra(i) = MAX(0., MIN(1. - dcf_con(i), cldfra(i) + dcf_mix(i)))974 qcld(i) = MAX(0., MIN(qtot(i) - dqi_con(i) - dqvc_con(i), qcld(i) + dqt_mix))841 cldfra(i) = MAX(0., MIN(1., cldfra(i) + dcf_mix(i))) 842 qcld(i) = MAX(0., MIN(qtot(i), qcld(i) + dqt_mix)) 975 843 qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dqvc_mix(i))) 976 844 977 ENDIF ! ( ( cldfra(i) .LT. ( 1. - dcf_con(i) - eps ) ) .AND. ( cldfra(i) .GT. eps ) ) 978 979 !--Finally, we add the tendencies of condensation 980 cldfra(i) = MIN(1., cldfra(i) + dcf_con(i)) 981 qcld(i) = MIN(qtot(i), qcld(i) + dqvc_con(i) + dqi_con(i)) 982 qvc(i) = MIN(qcld(i), qvc(i) + dqvc_con(i)) 845 ENDIF ! ( ( cldfra(i) .LT. ( 1. - eps ) ) .AND. ( cldfra(i) .GT. eps ) ) 983 846 984 847 … … 1052 915 !********************************************************************************** 1053 916 917 SUBROUTINE deposition_sublimation( & 918 qvapincld, qiceincld, temp, qsat, pplay, dtime, & 919 qvapincld_new) 920 921 USE lmdz_lscp_ini, ONLY: RV, RLSTT, RTT, RD, EPS_W 922 USE lmdz_lscp_ini, ONLY: eps 923 USE lmdz_lscp_ini, ONLY: depo_coef_cirrus, capa_cond_cirrus, rho_ice 924 925 REAL, INTENT(IN) :: qvapincld ! 926 REAL, INTENT(IN) :: qiceincld ! 927 REAL, INTENT(IN) :: temp ! 928 REAL, INTENT(IN) :: qsat ! 929 REAL, INTENT(IN) :: pplay ! 930 REAL, INTENT(IN) :: dtime ! time step [s] 931 REAL, INTENT(OUT):: qvapincld_new ! 932 933 ! 934 ! for unadjusted clouds 935 REAL :: qice_ratio 936 REAL :: pres_sat, rho, kappa 937 REAL :: air_thermal_conduct, water_vapor_diff 938 REAL :: iwc 939 REAL :: iwc_log_inf100, iwc_log_sup100 940 REAL :: iwc_inf100, alpha_inf100, coef_inf100 941 REAL :: mu_sup100, sigma_sup100, coef_sup100 942 REAL :: Dm_ice, rm_ice 943 944 !--If ok_unadjusted_clouds is set to TRUE, then the saturation adjustment 945 !--hypothesis is lost, and the vapor in the cloud is purely prognostic. 946 ! 947 !--The deposition equation is 948 !-- dmi/dt = alpha*4pi*C*Svi / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) ) 949 !--from Lohmann et al. (2016), where 950 !--alpha is the deposition coefficient [-] 951 !--mi is the mass of one ice crystal [kg] 952 !--C is the capacitance of an ice crystal [m] 953 !--Svi is the supersaturation ratio equal to (qvc - qsat)/qsat [-] 954 !--R_v is the specific gas constant for humid air [J/kg/K] 955 !--T is the temperature [K] 956 !--esi is the saturation pressure w.r.t. ice [Pa] 957 !--Dv is the diffusivity of water vapor [m2/s] 958 !--Ls is the specific latent heat of sublimation [J/kg/K] 959 !--ka is the thermal conductivity of dry air [J/m/s/K] 960 ! 961 !--alpha is a coefficient to take into account the fact that during deposition, a water 962 !--molecule cannot join the crystal from everywhere, it must do so that the crystal stays 963 !--coherent (with the same structure). It has no impact for sublimation. 964 !--We fix alpha = depo_coef_cirrus (=0.5 by default following Lohmann et al. (2016)) 965 !--during deposition, and alpha = 1. during sublimation. 966 !--The capacitance of the ice crystals is proportional to a parameter capa_cond_cirrus 967 !-- C = capa_cond_cirrus * rm_ice 968 ! 969 !--We have qice = Nice * mi, where Nice is the ice crystal 970 !--number concentration per kg of moist air 971 !--HYPOTHESIS 1: the ice crystals are spherical, therefore 972 !-- mi = 4/3 * pi * rm_ice**3 * rho_ice 973 !--HYPOTHESIS 2: the ice crystals are monodisperse with the 974 !--initial radius rm_ice_0. 975 !--NB. this is notably different than the assumption 976 !--of a distributed qice in the cloud made in the sublimation process; 977 !--should it be consistent? 978 ! 979 !--As the deposition process does not create new ice crystals, 980 !--and because we assume a same rm_ice value for all crystals 981 !--therefore the sublimation process does not destroy ice crystals 982 !--(or, in a limit case, it destroys all ice crystals), then 983 !--Nice is a constant during the sublimation/deposition process. 984 !-- dmi = dqi, et Nice = qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice ) 985 ! 986 !--The deposition equation then reads: 987 !-- dqi/dt = alpha*4pi*capa_cond_cirrus*rm_ice*(qvc-qsat)/qsat / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) ) * Nice 988 !-- dqi/dt = alpha*4pi*capa_cond_cirrus* (qi / qi_0)**(1/3) *rm_ice_0*(qvc-qsat)/qsat & 989 !-- / ( R_v*T/esi/Dv + Ls/ka/T * (Ls*R_v/T - 1) ) & 990 !-- * qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice ) 991 !-- dqi/dt = qi**(1/3) * (qvc - qsat) * qi_0**(2/3) & 992 !-- *alpha/qsat*capa_cond_cirrus/ (R_v*T/esi/Dv + Ls/ka/T*(Ls*R_v/T - 1)) / ( 1/3 rm_ice_0**2 rho_ice ) 993 !--and we have 994 !-- dqvc/dt = - qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2 995 !-- dqi/dt = qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2 996 !--where kappa = 1/3*rho_ice/capa_cond_cirrus*qsat*(R_v*T/esi/Dv + Ls/ka/T*(Ls/R_v/T - 1)) 997 ! 998 !--This system of equations can be resolved with an exact 999 !--explicit numerical integration, having one variable resolved 1000 !--explicitly, the other exactly. The exactly resolved variable is 1001 !--the one decreasing, so it is qvc if the process is 1002 !--condensation, qi if it is sublimation. 1003 ! 1004 !--kappa is computed as an initialisation constant, as it depends only 1005 !--on temperature and other pre-computed values 1006 pres_sat = qsat / ( EPS_W + ( 1. - EPS_W ) * qsat ) * pplay 1007 !--This formula for air thermal conductivity comes from Beard and Pruppacher (1971) 1008 air_thermal_conduct = ( 5.69 + 0.017 * ( temp - RTT ) ) * 1.e-3 * 4.184 1009 !--This formula for water vapor diffusivity comes from Hall and Pruppacher (1976) 1010 water_vapor_diff = 0.211 * ( temp / RTT )**1.94 * ( 101325. / pplay ) * 1.e-4 1011 kappa = 1. / 3. * rho_ice / capa_cond_cirrus * qsat & 1012 * ( RV * temp / water_vapor_diff / pres_sat & 1013 + RLSTT / air_thermal_conduct / temp * ( RLSTT / RV / temp - 1. ) ) 1014 !--NB. the greater kappa, the lower the efficiency of the deposition/sublimation process 1015 1016 !--Dry density [kg/m3] 1017 rho = pplay / temp / RD 1018 1019 !--Here, the initial vapor in the cloud is qvapincld, and we compute 1020 !--the new vapor qvapincld_new 1021 1022 !--rm_ice formula from McFarquhar and Heymsfield (1997) 1023 iwc = qiceincld * rho * 1e3 1024 iwc_inf100 = MIN(iwc, 0.252 * iwc**0.837) 1025 iwc_log_inf100 = LOG10( MAX(eps, iwc_inf100) ) 1026 iwc_log_sup100 = LOG10( MAX(eps, iwc - iwc_inf100) ) 1027 1028 alpha_inf100 = - 4.99E-3 - 0.0494 * iwc_log_inf100 1029 coef_inf100 = iwc_inf100 * alpha_inf100**3 / 120. 1030 1031 mu_sup100 = ( 5.2 + 0.0013 * ( temp - RTT ) ) & 1032 + ( 0.026 - 1.2E-3 * ( temp - RTT ) ) * iwc_log_sup100 1033 sigma_sup100 = ( 0.47 + 2.1E-3 * ( temp - RTT ) ) & 1034 + ( 0.018 - 2.1E-4 * ( temp - RTT ) ) * iwc_log_sup100 1035 coef_sup100 = ( iwc - iwc_inf100 ) / EXP( 3. * mu_sup100 + 4.5 * sigma_sup100**2 ) 1036 1037 Dm_ice = ( 2. / alpha_inf100 * coef_inf100 + EXP( mu_sup100 + 0.5 * sigma_sup100**2 ) & 1038 * coef_sup100 ) / ( coef_inf100 + coef_sup100 ) 1039 rm_ice = Dm_ice / 2. * 1.E-6 1040 1041 IF ( qvapincld .GE. qsat ) THEN 1042 !--If the cloud is initially supersaturated 1043 !--Exact explicit integration (qvc exact, qice explicit) 1044 qvapincld_new = qsat + ( qvapincld - qsat ) & 1045 * EXP( - depo_coef_cirrus * dtime * qiceincld / kappa / rm_ice**2 ) 1046 ELSE 1047 !--If the cloud is initially subsaturated 1048 !--Exact explicit integration (qice exact, qvc explicit) 1049 !--The barrier is set so that the resulting vapor in cloud 1050 !--cannot be greater than qsat 1051 !--qice_ratio is the ratio between the new ice content and 1052 !--the old one, it is comprised between 0 and 1 1053 qice_ratio = ( 1. - 2. / 3. / kappa / rm_ice**2 * dtime * ( qsat - qvapincld ) ) 1054 1055 IF ( qice_ratio .LT. 0. ) THEN 1056 !--The new vapor in cloud is set to 0 so that the 1057 !--sublimation process does not activate 1058 qvapincld_new = 0. 1059 ELSE 1060 !--Else, the sublimation process is activated with the 1061 !--diagnosed new cloud water vapor 1062 !--The new vapor in the cloud is increased with the 1063 !--sublimated ice 1064 qvapincld_new = qvapincld + qiceincld * ( 1. - qice_ratio**1.5 ) 1065 !--The new vapor in the cloud cannot be greater than qsat 1066 qvapincld_new = MIN(qvapincld_new, qsat) 1067 ENDIF ! qice_ratio .LT. 0. 1068 ENDIF ! qvapincld .GT. qsat 1069 1070 END SUBROUTINE deposition_sublimation 1071 1054 1072 END MODULE lmdz_lscp_condensation
Note: See TracChangeset
for help on using the changeset viewer.
