Changeset 5227 for LMDZ6/branches/Amaury_dev
- Timestamp:
- Sep 25, 2024, 10:39:24 AM (4 months ago)
- Location:
- LMDZ6/branches/Amaury_dev
- Files:
-
- 4 edited
-
. (modified) (1 prop)
-
libf/phylmd/lmdz_lscp.F90 (modified) (2 diffs)
-
libf/phylmd/lmdz_lscp_condensation.f90 (modified) (22 diffs)
-
libf/phylmd/lmdz_lscp_ini.F90 (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/branches/Amaury_dev
- Property svn:mergeinfo changed
/LMDZ6/trunk merged: 5210-5211
- Property svn:mergeinfo changed
-
LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_lscp.F90
r5226 r5227 118 118 USE lmdz_lscp_ini, ONLY: RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG 119 119 USE lmdz_lscp_ini, ONLY: ok_poprecip 120 USE lmdz_lscp_ini, ONLY: ok_ external_lognormal, ok_ice_supersat, ok_unadjusted_clouds, iflag_icefrac120 USE lmdz_lscp_ini, ONLY: ok_ice_supersat, ok_unadjusted_clouds, iflag_icefrac 121 121 USE lmdz_abort_physic, ONLY: abort_physic 122 122 IMPLICIT NONE … … 960 960 961 961 962 !--If .TRUE., calls an externalised version of the generalised 963 !--lognormal condensation scheme (Bony and Emanuel 2001) 964 !--Numerically, convergence is conserved with this option 965 !--The objective is to simplify LSCP 966 ELSE IF (ok_external_lognormal) THEN 962 ELSE 963 !--generalisedlognormal condensation scheme (Bony and Emanuel 2001) 964 967 965 968 966 CALL condensation_lognormal(& 969 klon, Tbef(:), zq(:), zqs(:), gammasat(:), ratqs(:, k), & 970 keepgoing(:), rneb(:, k), zqn(:), qvc(:)) 971 972 ELSE !--Generalised lognormal (Bony and Emanuel 2001) 973 974 DO i = 1, klon !todoan : check if loop in i is needed 975 976 IF (keepgoing(i)) THEN 977 978 zpdf_sig(i) = ratqs(i, k) * zq(i) 979 zpdf_k(i) = -sqrt(log(1. + (zpdf_sig(i) / zq(i))**2)) 980 zpdf_delta(i) = log(zq(i) / (gammasat(i) * zqs(i))) 981 zpdf_a(i) = zpdf_delta(i) / (zpdf_k(i) * sqrt(2.)) 982 zpdf_b(i) = zpdf_k(i) / (2. * sqrt(2.)) 983 zpdf_e1(i) = zpdf_a(i) - zpdf_b(i) 984 zpdf_e1(i) = sign(min(ABS(zpdf_e1(i)), 5.), zpdf_e1(i)) 985 zpdf_e1(i) = 1. - erf(zpdf_e1(i)) 986 zpdf_e2(i) = zpdf_a(i) + zpdf_b(i) 987 zpdf_e2(i) = sign(min(ABS(zpdf_e2(i)), 5.), zpdf_e2(i)) 988 zpdf_e2(i) = 1. - erf(zpdf_e2(i)) 989 990 IF (zpdf_e1(i)<1.e-10) THEN 991 rneb(i, k) = 0. 992 zqn(i) = zqs(i) 993 !--AB grid-mean vapor in the cloud - we assume saturation adjustment 994 qvc(i) = 0. 995 ELSE 996 rneb(i, k) = 0.5 * zpdf_e1(i) 997 zqn(i) = zq(i) * zpdf_e2(i) / zpdf_e1(i) 998 !--AB grid-mean vapor in the cloud - we assume saturation adjustment 999 qvc(i) = rneb(i, k) * zqs(i) 1000 END IF 1001 1002 END IF ! keepgoing 1003 END DO ! loop on klon 1004 1005 END IF ! .NOT. ok_ice_supersat 967 klon, Tbef, zq, zqs, gammasat, ratqs(:, k), & 968 keepgoing, rneb(:, k), zqn, qvc) 969 970 971 ENDIF ! .NOT. ok_ice_supersat 1006 972 1007 973 DO i = 1, klon -
LMDZ6/branches/Amaury_dev/libf/phylmd/lmdz_lscp_condensation.f90
r5224 r5227 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 … … 143 143 REAL, INTENT(IN) , DIMENSION(klon) :: ratio_qi_qtot ! specific ice water content to total specific water ratio [-] 144 144 REAL, INTENT(IN) , DIMENSION(klon) :: shear ! vertical shear [s-1] 145 REAL, INTENT(IN) , DIMENSION(klon) :: pbl_eps ! 146 REAL, INTENT(IN) , DIMENSION(klon) :: cell_area ! 145 REAL, INTENT(IN) , DIMENSION(klon) :: pbl_eps ! 146 REAL, INTENT(IN) , DIMENSION(klon) :: cell_area ! 147 147 REAL, INTENT(IN) , DIMENSION(klon) :: temp ! temperature [K] 148 148 REAL, INTENT(IN) , DIMENSION(klon) :: qtot ! total specific humidity (without precip) [kg/kg] … … 154 154 ! Input for aviation 155 155 ! 156 REAL, INTENT(IN), DIMENSION(klon) :: flight_dist ! 157 REAL, INTENT(IN), DIMENSION(klon) :: flight_h2o ! 156 REAL, INTENT(IN), DIMENSION(klon) :: flight_dist ! 157 REAL, INTENT(IN), DIMENSION(klon) :: flight_h2o ! 158 158 ! 159 159 ! Output … … 189 189 ! 190 190 REAL, INTENT(INOUT), DIMENSION(klon) :: Tcontr ! critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2) [K] 191 REAL, INTENT(INOUT), DIMENSION(klon) :: qcontr ! 192 REAL, INTENT(INOUT), DIMENSION(klon) :: qcontr2 ! 193 REAL, INTENT(INOUT), DIMENSION(klon) :: fcontrN ! 194 REAL, INTENT(INOUT), DIMENSION(klon) :: fcontrP ! 191 REAL, INTENT(INOUT), DIMENSION(klon) :: qcontr ! 192 REAL, INTENT(INOUT), DIMENSION(klon) :: qcontr2 ! 193 REAL, INTENT(INOUT), DIMENSION(klon) :: fcontrN ! 194 REAL, INTENT(INOUT), DIMENSION(klon) :: fcontrP ! 195 195 REAL, INTENT(INOUT), DIMENSION(klon) :: dcf_avi ! cloud fraction tendency because of aviation [s-1] 196 196 REAL, INTENT(INOUT), DIMENSION(klon) :: dqi_avi ! specific ice content tendency because of aviation [kg/kg/s] … … 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 … … 247 251 ! 248 252 !--more local variables for diagnostics 249 !--imported from YOMCST.h 253 !--imported from YOMCST.h 250 254 !--eps_w = 0.622 = ratio of molecular masses of water and dry air (kg H2O kg air -1) 251 255 !--RCPD = 1004 J kg air−1 K−1 = the isobaric heat capacity of air … … 254 258 !--Qheat = 43. / 50. / 120. MJ / kg fuel for kerosene / methane / dihydrogen 255 259 !REAL, PARAMETER :: EiH2O=1.25 !--emission index of water vapour for kerosene (kg kg-1) 256 !REAL, PARAMETER :: Qheat=43.E6 !--specific combustion heat for kerosene (J kg-1) 260 !REAL, PARAMETER :: Qheat=43.E6 !--specific combustion heat for kerosene (J kg-1) 257 261 !REAL, PARAMETER :: eta=0.3 !--average propulsion efficiency of the aircraft 258 !--Gcontr is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute 262 !--Gcontr is the slope of the mean phase trajectory in the turbulent exhaust field on an absolute 259 263 !--temperature versus water vapor partial pressure diagram. G has the unit of Pa K−1. Rap et al JGR 2010. 260 264 !REAL :: Gcontr 261 !--Tcontr = critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2) 265 !--Tcontr = critical temperature for contrail formation (T_LM in Schumann 1996, Eq 31 in appendix 2) 262 266 !--qsatliqcontr = e_L(T_LM) in Schumann 1996 but expressed in specific humidity (kg kg humid air-1) 263 267 !REAL :: qsatliqcontr … … 270 274 ! Ajout des émissions de H2O dues à l'aviation 271 275 ! q is the specific humidity (kg/kg humid air) hence the complicated equation to update q 272 ! qnew = ( m_humid_air * qold + dm_H2O ) / ( m_humid_air + dm_H2O ) 273 ! = ( m_dry_air * qold + dm_h2O * (1-qold) ) / (m_dry_air + dm_H2O * (1-qold) ) 274 ! The equation is derived by writing m_humid_air = m_dry_air + m_H2O = m_dry_air / (1-q) 276 ! qnew = ( m_humid_air * qold + dm_H2O ) / ( m_humid_air + dm_H2O ) 277 ! = ( m_dry_air * qold + dm_h2O * (1-qold) ) / (m_dry_air + dm_H2O * (1-qold) ) 278 ! The equation is derived by writing m_humid_air = m_dry_air + m_H2O = m_dry_air / (1-q) 275 279 ! flight_h2O is in kg H2O / s / cell 276 280 ! 277 !IF (ok_plane_h2o) THEN 278 ! q = ( M_cell*q + flight_h2o(i,k)*dtime*(1.-q) ) / (M_cell + flight_h2o(i,k)*dtime*(1.-q) ) 281 !IF (ok_plane_h2o) THEN 282 ! q = ( M_cell*q + flight_h2o(i,k)*dtime*(1.-q) ) / (M_cell + flight_h2o(i,k)*dtime*(1.-q) ) 279 283 !ENDIF 280 284 … … 345 349 !--can be greater than the total water in the gridbox 346 350 cldfra(i) = MAX(0., MIN(1., cf_seri(i))) 347 qcld(i) = MAX(0., MIN(qtot(i), ( ratio_qi_qtot(i) + rvc_seri(i) ) * qtot(i))) 351 qcld(i) = MAX(0., MIN(qtot(i), ( ratio_qi_qtot(i) + rvc_seri(i) ) * qtot(i))) 348 352 qvc(i) = MAX(0., MIN(qcld(i), rvc_seri(i) * qtot(i))) 349 353 ok_warm_cloud = .FALSE. … … 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 ! … … 458 462 !-- SUBLIMATION OF ICE AND DEPOSITION OF VAPOR IN THE CLOUD -- 459 463 !------------------------------------------------------------------- 460 464 461 465 !--If there is a cloud 462 466 IF ( cldfra(i) > eps ) THEN 463 467 464 468 qvapincld = qvc(i) / cldfra(i) 465 469 qiceincld = ( qcld(i) / cldfra(i) - qvapincld ) 466 470 467 471 !--If the ice water content is too low, the cloud is purely sublimated 468 472 !--Most probably, we advected a cloud with no ice water content (possible … … 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 >= 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 < 0. ) THEN … … 539 558 ENDIF ! ok_unadjusted_clouds 540 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 !------------------------------------ 574 541 575 !--If the vapor in cloud is below vapor needed for the cloud to survive 542 576 IF ( qvapincld < qvapincld_new ) THEN … … 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 > 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) < rhlmid_pdf_lscp ) .OR. ( pdf_e4 > 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 682 685 ELSEIF ( cf_cond > eps ) THEN 686 687 dcf_con(i) = ( 1. - cldfra(i) ) * cf_cond 688 dqt_con = ( 1. - cldfra(i) ) * qt_cond 689 683 690 !--Barriers 684 691 dcf_con(i) = MIN(dcf_con(i), 1. - cldfra(i)) … … 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/branches/Amaury_dev/libf/phylmd/lmdz_lscp_ini.F90
r5224 r5227 144 144 145 145 !--Parameters for condensation and ice supersaturation 146 LOGICAL, SAVE, PROTECTED :: ok_external_lognormal=.FALSE. ! if True, the lognormal condensation scheme is calculated in the lmdz_lscp_condensation routine147 !$OMP THREADPRIVATE(ok_external_lognormal)148 146 149 147 LOGICAL, SAVE, PROTECTED :: ok_ice_supersat=.FALSE. ! activates the condensation scheme that allows for ice supersaturation … … 186 184 !$OMP THREADPRIVATE(rhl0_pdf_lscp) 187 185 188 REAL, SAVE, PROTECTED :: cond_thresh_pdf_lscp=1.E-10 ! [-] threshold for the formation of new cloud 189 !$OMP THREADPRIVATE(cond_thresh_pdf_lscp) 190 186 191 187 REAL, SAVE, PROTECTED :: a_homofreez=2.349 ! [-] factor for the Koop homogeneous freezing fit 192 188 !$OMP THREADPRIVATE(a_homofreez) … … 421 417 CALL getin_p('snow_fallspeed_cld',snow_fallspeed_cld) 422 418 ! for condensation and ice supersaturation 423 CALL getin_p('ok_external_lognormal',ok_external_lognormal)424 419 CALL getin_p('ok_unadjusted_clouds',ok_unadjusted_clouds) 425 420 CALL getin_p('ok_weibull_warm_clouds',ok_weibull_warm_clouds) … … 434 429 CALL getin_p('kappa_pdf_lscp',kappa_pdf_lscp) 435 430 CALL getin_p('rhl0_pdf_lscp',rhl0_pdf_lscp) 436 CALL getin_p('cond_thresh_pdf_lscp',cond_thresh_pdf_lscp)437 431 CALL getin_p('a_homofreez',a_homofreez) 438 432 CALL getin_p('b_homofreez',b_homofreez) … … 503 497 WRITE(lunout,*) 'lscp_ini, snow_fallspeed_cld:', snow_fallspeed_cld 504 498 ! for condensation and ice supersaturation 505 WRITE(lunout,*) 'lscp_ini, ok_external_lognormal:', ok_external_lognormal506 499 WRITE(lunout,*) 'lscp_ini, ok_ice_supersat:', ok_ice_supersat 507 500 WRITE(lunout,*) 'lscp_ini, ok_unadjusted_clouds:', ok_unadjusted_clouds … … 543 536 544 537 545 !--Check flags for condensation and ice supersaturation546 IF ( ok_external_lognormal .AND. ok_ice_supersat ) THEN547 abort_message = 'in lscp, ok_external_lognormal=y is incompatible with ok_ice_supersat=y'548 CALL abort_physic (modname,abort_message,1)549 ENDIF550 551 538 IF ( ok_weibull_warm_clouds .AND. .NOT. ok_ice_supersat ) THEN 552 539 abort_message = 'in lscp, ok_weibull_warm_clouds=y needs ok_ice_supersat=y'
Note: See TracChangeset
for help on using the changeset viewer.
