- Timestamp:
- Sep 21, 2024, 11:17:30 AM (5 hours ago)
- Location:
- LMDZ6/trunk/libf/phylmd
- Files:
-
- 3 edited
-
lmdz_lscp.F90 (modified) (1 diff)
-
lmdz_lscp_condensation.F90 (modified) (14 diffs)
-
lmdz_lscp_ini.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/lmdz_lscp.F90
r5208 r5210 975 975 976 976 CALL condensation_lognormal( & 977 klon, Tbef (:), zq(:), zqs(:), gammasat(:), ratqs(:,k), &978 keepgoing (:), rneb(:,k), zqn(:), qvc(:))977 klon, Tbef, zq, zqs, gammasat, ratqs(:,k), & 978 keepgoing, rneb(:,k), zqn, qvc) 979 979 980 980 -
LMDZ6/trunk/libf/phylmd/lmdz_lscp_condensation.F90
r5204 r5210 124 124 mu_subl_pdf_lscp, beta_pdf_lscp, temp_thresh_pdf_lscp, & 125 125 rhlmid_pdf_lscp, k0_pdf_lscp, kappa_pdf_lscp, rhl0_pdf_lscp, & 126 co nd_thresh_pdf_lscp, coef_mixing_lscp, coef_shear_lscp, &126 coef_mixing_lscp, coef_shear_lscp, & 127 127 chi_mixing_lscp, rho_ice 128 128 … … 212 212 REAL :: pres_sat, kappa 213 213 REAL :: air_thermal_conduct, water_vapor_diff 214 REAL :: r_ice 214 REAL :: iwc 215 REAL :: iwc_log_inf100, iwc_log_sup100 216 REAL :: iwc_inf100, alpha_inf100, coef_inf100 217 REAL :: mu_sup100, sigma_sup100, coef_sup100 218 REAL :: Dm_ice, rm_ice 215 219 ! 216 220 ! for sublimation … … 404 408 !--during deposition, and alpha = 1. during sublimation. 405 409 !--The capacitance of the ice crystals is proportional to a parameter capa_cond_cirrus 406 !-- C = capa_cond_cirrus * r _ice410 !-- C = capa_cond_cirrus * rm_ice 407 411 ! 408 412 !--We have qice = Nice * mi, where Nice is the ice crystal 409 413 !--number concentration per kg of moist air 410 414 !--HYPOTHESIS 1: the ice crystals are spherical, therefore 411 !-- mi = 4/3 * pi * r _ice**3 * rho_ice415 !-- mi = 4/3 * pi * rm_ice**3 * rho_ice 412 416 !--HYPOTHESIS 2: the ice crystals are monodisperse with the 413 !--initial radius r _ice_0.417 !--initial radius rm_ice_0. 414 418 !--NB. this is notably different than the assumption 415 419 !--of a distributed qice in the cloud made in the sublimation process; … … 417 421 ! 418 422 !--As the deposition process does not create new ice crystals, 419 !--and because we assume a same r _ice value for all crystals423 !--and because we assume a same rm_ice value for all crystals 420 424 !--therefore the sublimation process does not destroy ice crystals 421 425 !--(or, in a limit case, it destroys all ice crystals), then 422 426 !--Nice is a constant during the sublimation/deposition process. 423 !-- dmi = dqi, et Nice = qi_0 / ( 4/3 RPI r _ice_0**3 rho_ice )427 !-- dmi = dqi, et Nice = qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice ) 424 428 ! 425 429 !--The deposition equation then reads: 426 !-- dqi/dt = alpha*4pi*capa_cond_cirrus*r _ice*(qvc-qsat)/qsat / ( R_v*T/esi/Dv + Ls/ka/T * (Ls/R_v/T - 1) ) * Nice427 !-- dqi/dt = alpha*4pi*capa_cond_cirrus* (qi / qi_0)**(1/3) *r _ice_0*(qvc-qsat)/qsat &430 !-- 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 431 !-- dqi/dt = alpha*4pi*capa_cond_cirrus* (qi / qi_0)**(1/3) *rm_ice_0*(qvc-qsat)/qsat & 428 432 !-- / ( R_v*T/esi/Dv + Ls/ka/T * (Ls*R_v/T - 1) ) & 429 !-- * qi_0 / ( 4/3 RPI r _ice_0**3 rho_ice )433 !-- * qi_0 / ( 4/3 RPI rm_ice_0**3 rho_ice ) 430 434 !-- dqi/dt = qi**(1/3) * (qvc - qsat) * qi_0**(2/3) & 431 !-- *alpha/qsat*capa_cond_cirrus/ (R_v*T/esi/Dv + Ls/ka/T*(Ls*R_v/T - 1)) / ( 1/3 r _ice_0**2 rho_ice )435 !-- *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 ) 432 436 !--and we have 433 !-- dqvc/dt = - qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / r _ice_0**2434 !-- dqi/dt = qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / r _ice_0**2437 !-- dqvc/dt = - qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2 438 !-- dqi/dt = qi**(1/3) * (qvc - qsat) / kappa * alpha * qi_0**(2/3) / rm_ice_0**2 435 439 !--where kappa = 1/3*rho_ice/capa_cond_cirrus*qsat*(R_v*T/esi/Dv + Ls/ka/T*(Ls/R_v/T - 1)) 436 440 ! … … 489 493 !--the new vapor qvapincld_new 490 494 491 !--r_ice formula from Sun and Rikus (1999) 492 r_ice = 1.e-6 * ( 45.8966 * ( qiceincld * rho * 1e3 )**0.2214 & 493 + 0.7957 * ( qiceincld * rho * 1e3 )**0.2535 * ( temp(i) - RTT + 190. ) ) / 2. 495 !--rm_ice formula from McFarquhar and Heymsfield (1997) 496 iwc = qiceincld * rho * 1e3 497 iwc_inf100 = MIN(iwc, 0.252 * iwc**0.837) 498 iwc_log_inf100 = LOG10( MAX(eps, iwc_inf100) ) 499 iwc_log_sup100 = LOG10( MAX(eps, iwc - iwc_inf100) ) 500 501 alpha_inf100 = - 4.99E-3 - 0.0494 * iwc_log_inf100 502 coef_inf100 = iwc_inf100 * alpha_inf100**3. / 120. 503 504 mu_sup100 = ( 5.2 + 0.0013 * ( temp(i) - RTT ) ) & 505 + ( 0.026 - 1.2E-3 * ( temp(i) - RTT ) ) * iwc_log_sup100 506 sigma_sup100 = ( 0.47 + 2.1E-3 * ( temp(i) - RTT ) ) & 507 + ( 0.018 - 2.1E-4 * ( temp(i) - RTT ) ) * iwc_log_sup100 508 coef_sup100 = ( iwc - iwc_inf100 ) / EXP( 3 * mu_sup100 + 4.5 * sigma_sup100**2. ) 509 510 Dm_ice = ( 2. / alpha_inf100 * coef_inf100 + EXP( mu_sup100 + 0.5 * sigma_sup100**2. ) & 511 * coef_sup100 ) / ( coef_inf100 + coef_sup100 ) 512 rm_ice = Dm_ice / 2. * 1.E-6 494 513 495 514 IF ( qvapincld .GE. qsat(i) ) THEN … … 497 516 !--Exact explicit integration (qvc exact, qice explicit) 498 517 qvapincld_new = qsat(i) + ( qvapincld - qsat(i) ) & 499 * EXP( - depo_coef_cirrus * dtime * qiceincld / kappa / r _ice**2. )518 * EXP( - depo_coef_cirrus * dtime * qiceincld / kappa / rm_ice**2. ) 500 519 ELSE 501 520 !--If the cloud is initially subsaturated … … 505 524 !--qice_ratio is the ratio between the new ice content and 506 525 !--the old one, it is comprised between 0 and 1 507 qice_ratio = ( 1. - 2. / 3. / kappa / r _ice**2. * dtime * ( qsat(i) - qvapincld ) )526 qice_ratio = ( 1. - 2. / 3. / kappa / rm_ice**2. * dtime * ( qsat(i) - qvapincld ) ) 508 527 509 528 IF ( qice_ratio .LT. 0. ) THEN … … 538 557 qvapincld_new = qsat(i) 539 558 ENDIF ! ok_unadjusted_clouds 559 560 !--Adjustment of the IWC to the new vapor in cloud 561 !--(this can be either positive or negative) 562 dqvc_adj(i) = ( qvapincld_new * cldfra(i) - qvc(i) ) 563 dqi_adj(i) = - dqvc_adj(i) 564 565 !--Add tendencies 566 !--The vapor in the cloud is updated, but not qcld as it is constant 567 !--through this process, as well as cldfra which is unmodified 568 qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dqvc_adj(i))) 569 570 571 !------------------------------------ 572 !-- DISSIPATION OF THE CLOUD -- 573 !------------------------------------ 540 574 541 575 !--If the vapor in cloud is below vapor needed for the cloud to survive … … 577 611 578 612 !--Tendencies and diagnostics 579 dqvc_sub(i) = dcf_sub(i) * qvapincld 580 dqi_sub(i) = dqt_sub - dqvc_sub(i) 613 dqvc_sub(i) = dqt_sub 581 614 582 615 !--Add tendencies … … 586 619 587 620 ENDIF ! qvapincld .LT. qvapincld_new 588 589 !--Adjustment of the IWC of the remaining cloud590 !--to the new vapor in cloud (this can be either positive or negative)591 dqvc_adj(i) = ( qvapincld_new * cldfra(i) - qvc(i) )592 dqi_adj(i) = - dqvc_adj(i)593 594 !--Add tendencies595 !--The vapor in the cloud is updated, but not qcld as it is constant596 !--through this process, as well as cldfra which is unmodified597 qvc(i) = MAX(0., MIN(qcld(i), qvc(i) + dqvc_adj(i)))598 621 599 622 ENDIF ! qiceincld .LT. eps … … 653 676 654 677 IF ( ok_warm_cloud ) THEN 655 !--If the statistical scheme is activated, the first656 !-- calculation of the tendencies is kept for the diagnosis of657 !-- the cloud properties, and thecondensation process stops678 !--If the statistical scheme is activated, the calculated increase is equal 679 !--to the cloud fraction, we assume saturation adjustment, and the 680 !--condensation process stops 658 681 cldfra(i) = cf_cond 659 682 qcld(i) = qt_cond 660 683 qvc(i) = cldfra(i) * qsat(i) 661 684 662 ELSEIF ( cf_cond .GT. cond_thresh_pdf_lscp ) THEN 663 !--We check that there is something to condense, if not the 664 !--condensation process stops 665 666 !--Calculation of the limit qclr in cloud value (stored in pdf_e4) 667 668 pdf_e3 = ( LOG( 1. / cond_thresh_pdf_lscp ) ** ( 1. / pdf_alpha ) - pdf_e2 ) & 669 / SQRT( GAMMA(1. + 2. / pdf_alpha) - pdf_e2**2. ) 670 pdf_e4 = ( 2. * pdf_e3 * pdf_e1 - pdf_x ) & 671 / ( pdf_e3 * pdf_e1 / rhlmid_pdf_lscp - 1. ) 672 IF ( ( MIN(pdf_e4, rhl_clr) .LT. rhlmid_pdf_lscp ) .OR. ( pdf_e4 .GT. rhl_clr ) ) THEN 673 pdf_e4 = pdf_x / ( pdf_e3 * pdf_e1 / rhlmid_pdf_lscp + 1. ) 674 ENDIF 675 pdf_e4 = pdf_e4 * qsatl(i) / 100. 676 677 !--Calculation of the final tendencies 678 dcf_con(i) = ( 1. - cldfra(i) ) * ( pdf_e4 - qclr / ( 1. - cldfra(i) ) ) & 679 / ( pdf_e4 - qt_cond / cf_cond ) 680 dqt_con = qclr - pdf_e4 * ( 1. - cldfra(i) - dcf_con(i) ) 681 685 ELSEIF ( cf_cond .GT. eps ) THEN 686 687 dcf_con(i) = ( 1. - cldfra(i) ) * cf_cond 688 dqt_con = ( 1. - cldfra(i) ) * qt_cond 682 689 683 690 !--Barriers … … 692 699 qvapincld = gamma_cond(i) * qsat(i) 693 700 qiceincld = dqt_con / dcf_con(i) - gamma_cond(i) * qsat(i) 694 !--r_ice formula from Sun and Rikus (1999) 695 r_ice = 1.e-6 * ( 45.8966 * ( qiceincld * rho * 1e3 )**0.2214 & 696 + 0.7957 * ( qiceincld * rho * 1e3 )**0.2535 * ( temp(i) - RTT + 190. ) ) / 2. 701 702 !--rm_ice formula from McFarquhar and Heymsfield (1997) 703 iwc = qiceincld * rho * 1e3 704 iwc_inf100 = MIN(iwc, 0.252 * iwc**0.837) 705 iwc_log_inf100 = LOG10( MAX(eps, iwc_inf100) ) 706 iwc_log_sup100 = LOG10( MAX(eps, iwc - iwc_inf100) ) 707 708 alpha_inf100 = - 4.99E-3 - 0.0494 * iwc_log_inf100 709 coef_inf100 = iwc_inf100 * alpha_inf100**3. / 120. 710 711 mu_sup100 = ( 5.2 + 0.0013 * ( temp(i) - RTT ) ) & 712 + ( 0.026 - 1.2E-3 * ( temp(i) - RTT ) ) * iwc_log_sup100 713 sigma_sup100 = ( 0.47 + 2.1E-3 * ( temp(i) - RTT ) ) & 714 + ( 0.018 - 2.1E-4 * ( temp(i) - RTT ) ) * iwc_log_sup100 715 coef_sup100 = ( iwc - iwc_inf100 ) / EXP( 3. * mu_sup100 + 4.5 * sigma_sup100**2. ) 716 717 Dm_ice = ( 2. / alpha_inf100 * coef_inf100 + EXP( mu_sup100 + 0.5 * sigma_sup100**2. ) & 718 * coef_sup100 ) / ( coef_inf100 + coef_sup100 ) 719 rm_ice = Dm_ice / 2. * 1.E-6 697 720 !--As qvapincld is necessarily greater than qsat, we only 698 721 !--use the exact explicit formulation 699 722 !--Exact explicit version 700 723 qvapincld = qsat(i) + ( qvapincld - qsat(i) ) & 701 * EXP( - depo_coef_cirrus * dtime / 2. * qiceincld / kappa / r _ice**2. )724 * EXP( - depo_coef_cirrus * dtime / 2. * qiceincld / kappa / rm_ice**2. ) 702 725 ELSE 703 726 !--We keep the saturation adjustment hypothesis, and the vapor in the … … 715 738 !-- qtot(i) - dqi_con(i) - dqvc_con(i) 716 739 717 ENDIF ! ok_warm_cloud, cf_cond .GT. cond_thresh_pdf_lscp740 ENDIF ! ok_warm_cloud, cf_cond .GT. eps 718 741 ENDIF ! ( 1. - cldfra(i) ) .GT. eps 719 742 … … 854 877 L_shear = coef_shear_lscp * shear(i) * dz * dtime 855 878 !--therefore, the total increase in fraction is 856 shear_fra = RPI * L_shear * a_mix * bovera * N_cld_mix & 879 !-- N_cld_mix * ( (a + L_shear) * b - a * b ) * RPI / 2. / cell_area 880 shear_fra = RPI * L_shear * a_mix * bovera / 2. * N_cld_mix & 857 881 / cell_area(i) / ( 1. - dcf_con(i) ) 858 882 !--and the environment and cloud mixed fractions are the same, -
LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.F90
r5204 r5210 185 185 REAL, SAVE, PROTECTED :: rhl0_pdf_lscp=88.7 ! [%] factor for the PDF fit of water vapor in UTLS 186 186 !$OMP THREADPRIVATE(rhl0_pdf_lscp) 187 188 REAL, SAVE, PROTECTED :: cond_thresh_pdf_lscp=1.E-10 ! [-] threshold for the formation of new cloud189 !$OMP THREADPRIVATE(cond_thresh_pdf_lscp)190 187 191 188 REAL, SAVE, PROTECTED :: a_homofreez=2.349 ! [-] factor for the Koop homogeneous freezing fit … … 433 430 CALL getin_p('kappa_pdf_lscp',kappa_pdf_lscp) 434 431 CALL getin_p('rhl0_pdf_lscp',rhl0_pdf_lscp) 435 CALL getin_p('cond_thresh_pdf_lscp',cond_thresh_pdf_lscp)436 432 CALL getin_p('a_homofreez',a_homofreez) 437 433 CALL getin_p('b_homofreez',b_homofreez)
Note: See TracChangeset
for help on using the changeset viewer.
