Changeset 5996
- Timestamp:
- Jan 5, 2026, 5:27:37 PM (7 days ago)
- Location:
- LMDZ6/trunk/libf
- Files:
-
- 2 added
- 4 edited
-
phylmd/lmdz_call_cloud_optics_prop.f90 (modified) (1 diff)
-
phylmd/lmdz_lscp_main.f90 (modified) (1 diff)
-
phylmd/lmdz_lscp_phase.f90 (added)
-
phylmd/lmdz_lscp_tools.f90 (modified) (1 diff)
-
phylmd/nuage.f90 (modified) (1 diff)
-
phylmdiso/lmdz_lscp_phase.f90 (added)
Legend:
- Unmodified
- Added
- Removed
-
LMDZ6/trunk/libf/phylmd/lmdz_call_cloud_optics_prop.f90
r5828 r5996 39 39 40 40 USE icefrac_lsc_mod ! computes ice fraction (JBM 3/14) 41 USE lmdz_lscp_ tools, only: icefrac_lscp41 USE lmdz_lscp_phase, only: icefrac_lscp 42 42 USE nuage_mod, ONLY: nuage 43 43 -
LMDZ6/trunk/libf/phylmd/lmdz_lscp_main.f90
r5994 r5996 71 71 72 72 ! USE de modules contenant des fonctions. 73 USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, calc_gammasat 74 USE lmdz_lscp_ tools, ONLY : icefrac_lscp, icefrac_lscp_turb, distance_to_cloud_top73 USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, calc_gammasat, distance_to_cloud_top 74 USE lmdz_lscp_phase, ONLY : icefrac_lscp, icefrac_lscp_turb 75 75 USE lmdz_lscp_condensation, ONLY : condensation_lognormal, condensation_ice_supersat 76 76 USE lmdz_lscp_condensation, ONLY : cloudth, cloudth_v3, cloudth_v6, condensation_cloudth -
LMDZ6/trunk/libf/phylmd/lmdz_lscp_tools.f90
r5933 r5996 100 100 END SUBROUTINE FALLICE_VELOCITY 101 101 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 102 103 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++104 SUBROUTINE ICEFRAC_LSCP(klon, temp, iflag_ice_thermo, distcltop, temp_cltop, icefrac, dicefracdT)105 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++106 107 ! Compute the ice fraction 1-xliq (see e.g.108 ! Doutriaux-Boucher & Quaas 2004, section 2.2.)109 ! as a function of temperature110 ! see also Fig 3 of Madeleine et al. 2020, JAMES111 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++112 113 114 USE print_control_mod, ONLY: lunout, prt_level115 USE lmdz_lscp_ini, ONLY: t_glace_min, t_glace_max, exposant_glace, iflag_t_glace116 USE lmdz_lscp_ini, ONLY : RTT, dist_liq, temp_nowater117 118 IMPLICIT NONE119 120 121 INTEGER, INTENT(IN) :: klon ! number of horizontal grid points122 REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature123 REAL, INTENT(IN), DIMENSION(klon) :: distcltop ! distance to cloud top124 REAL, INTENT(IN), DIMENSION(klon) :: temp_cltop ! temperature of cloud top125 INTEGER, INTENT(IN) :: iflag_ice_thermo126 REAL, INTENT(OUT), DIMENSION(klon) :: icefrac127 REAL, INTENT(OUT), DIMENSION(klon) :: dicefracdT128 129 130 INTEGER i131 REAL liqfrac_tmp, dicefrac_tmp132 REAL Dv, denomdep,beta,qsi,dqsidt133 LOGICAL ice_thermo134 135 CHARACTER (len = 20) :: modname = 'lscp_tools'136 CHARACTER (len = 80) :: abort_message137 138 IF ((iflag_t_glace.LT.2)) THEN !.OR. (iflag_t_glace.GT.6)) THEN139 abort_message = 'lscp cannot be used if iflag_t_glace<2 or >6'140 CALL abort_physic(modname,abort_message,1)141 ENDIF142 143 IF (.NOT.((iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3))) THEN144 abort_message = 'lscp cannot be used without ice thermodynamics'145 CALL abort_physic(modname,abort_message,1)146 ENDIF147 148 149 DO i=1,klon150 151 ! old function with sole dependence upon temperature152 IF (iflag_t_glace .EQ. 2) THEN153 liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)154 liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)155 icefrac(i) = (1.0-liqfrac_tmp)**exposant_glace156 IF (icefrac(i) .GT.0.) THEN157 dicefracdT(i)= exposant_glace * (icefrac(i)**(exposant_glace-1.)) &158 / (t_glace_min - t_glace_max)159 ENDIF160 161 IF ((icefrac(i).EQ.0).OR.(icefrac(i).EQ.1)) THEN162 dicefracdT(i)=0.163 ENDIF164 165 ENDIF166 167 ! function of temperature used in CMIP6 physics168 IF (iflag_t_glace .EQ. 3) THEN169 liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)170 liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)171 icefrac(i) = 1.0-liqfrac_tmp**exposant_glace172 IF ((icefrac(i) .GT.0.) .AND. (liqfrac_tmp .GT. 0.)) THEN173 dicefracdT(i)= exposant_glace * ((liqfrac_tmp)**(exposant_glace-1.)) &174 / (t_glace_min - t_glace_max)175 ELSE176 dicefracdT(i)=0.177 ENDIF178 ENDIF179 180 ! for iflag_t_glace .GE. 4, the liquid fraction depends upon temperature at cloud top181 ! and then decreases with decreasing height182 183 !with linear function of temperature at cloud top184 IF (iflag_t_glace .EQ. 4) THEN185 liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)186 liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)187 icefrac(i) = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)188 dicefrac_tmp = - temp(i)/(t_glace_max-t_glace_min)189 dicefracdT(i) = dicefrac_tmp*exp(-distcltop(i)/dist_liq)190 IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN191 dicefracdT(i) = 0.192 ENDIF193 ENDIF194 195 ! with CMIP6 function of temperature at cloud top196 IF ((iflag_t_glace .EQ. 5) .OR. (iflag_t_glace .EQ. 7)) THEN197 liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min)198 liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)199 liqfrac_tmp = liqfrac_tmp**exposant_glace200 icefrac(i) = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)201 IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN202 dicefracdT(i) = 0.203 ELSE204 dicefracdT(i) = exposant_glace*((liqfrac_tmp)**(exposant_glace-1.))/(t_glace_min- t_glace_max) &205 *exp(-distcltop(i)/dist_liq)206 ENDIF207 ENDIF208 209 ! with modified function of temperature at cloud top210 ! to get largere values around 260 K, works well with t_glace_min = 241K211 IF (iflag_t_glace .EQ. 6) THEN212 IF (temp(i) .GT. t_glace_max) THEN213 liqfrac_tmp = 1.214 ELSE215 liqfrac_tmp = -((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))**2+1.216 ENDIF217 liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0)218 icefrac(i) = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.)219 IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN220 dicefracdT(i) = 0.221 ELSE222 dicefracdT(i) = 2*((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))/(t_glace_max-t_glace_min) &223 *exp(-distcltop(i)/dist_liq)224 ENDIF225 ENDIF226 227 ! if temperature or temperature of cloud top <-40°C,228 IF (iflag_t_glace .GE. 4) THEN229 IF ((temp_cltop(i) .LE. temp_nowater) .AND. (temp(i) .LE. t_glace_max)) THEN230 icefrac(i) = 1.231 dicefracdT(i) = 0.232 ENDIF233 ENDIF234 235 236 ENDDO ! klon237 RETURN238 239 END SUBROUTINE ICEFRAC_LSCP240 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++241 242 243 SUBROUTINE ICEFRAC_LSCP_TURB(klon, dtime, pticefracturb, temp, pplay, paprsdn, paprsup, wvel, qice_ini, snowcld, qtot_incl, cldfra, tke, &244 tke_dissip, sursat_e, invtau_e, qliq, qvap_cld, qice, icefrac, dicefracdT, cldfraliq, sigma2_icefracturb, mean_icefracturb)245 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++246 ! Compute the liquid, ice and vapour content (+ice fraction) based247 ! on turbulence (see Fields 2014, Furtado 2016, Raillard 2025)248 ! L.Raillard (23/09/24)249 ! E.Vignon (03/2025) : additional elements for treatment of convective250 ! boundary layer clouds251 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++252 253 254 USE lmdz_lscp_ini, ONLY : prt_level, lunout255 USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI256 USE lmdz_lscp_ini, ONLY : seuil_neb, temp_nowater257 USE lmdz_lscp_ini, ONLY : naero5, gamma_snwretro, gamma_taud, capa_crystal, rho_ice258 USE lmdz_lscp_ini, ONLY : eps, snow_fallspeed259 260 IMPLICIT NONE261 262 INTEGER, INTENT(IN) :: klon !--number of horizontal grid points263 REAL, INTENT(IN) :: dtime !--time step [s]264 LOGICAL, INTENT(IN), DIMENSION(klon) :: pticefracturb !--grid points concerned by this routine265 REAL, INTENT(IN), DIMENSION(klon) :: temp !--temperature266 REAL, INTENT(IN), DIMENSION(klon) :: pplay !--pressure in the middle of the layer [Pa]267 REAL, INTENT(IN), DIMENSION(klon) :: paprsdn !--pressure at the bottom interface of the layer [Pa]268 REAL, INTENT(IN), DIMENSION(klon) :: paprsup !--pressure at the top interface of the layer [Pa]269 REAL, INTENT(IN), DIMENSION(klon) :: wvel !--vertical velocity [m/s]270 REAL, INTENT(IN), DIMENSION(klon) :: qtot_incl !--specific total cloud water in-cloud content [kg/kg]271 REAL, INTENT(IN), DIMENSION(klon) :: cldfra !--cloud fraction in gridbox [-]272 REAL, INTENT(IN), DIMENSION(klon) :: tke !--turbulent kinetic energy [m2/s2]273 REAL, INTENT(IN), DIMENSION(klon) :: tke_dissip !--TKE dissipation [m2/s3]274 275 REAL, INTENT(IN), DIMENSION(klon) :: qice_ini !--initial specific ice content gridbox-mean [kg/kg]276 REAL, INTENT(IN), DIMENSION(klon) :: snowcld !--in-cloud snowfall flux [kg/m2/s]277 REAL, INTENT(IN), DIMENSION(klon) :: sursat_e !--environment supersaturation [-]278 REAL, INTENT(IN), DIMENSION(klon) :: invtau_e !--inverse time-scale of mixing with environment [s-1]279 REAL, INTENT(OUT), DIMENSION(klon) :: qliq !--specific liquid content gridbox-mean [kg/kg]280 REAL, INTENT(OUT), DIMENSION(klon) :: qvap_cld !--specific cloud vapor content, gridbox-mean [kg/kg]281 REAL, INTENT(OUT), DIMENSION(klon) :: qice !--specific ice content gridbox-mean [kg/kg]282 283 REAL, INTENT(INOUT), DIMENSION(klon) :: icefrac !--fraction of ice in condensed water [-]284 REAL, INTENT(INOUT), DIMENSION(klon) :: dicefracdT285 286 REAL, INTENT(OUT), DIMENSION(klon) :: cldfraliq !--fraction of cldfra where liquid [-]287 REAL, INTENT(OUT), DIMENSION(klon) :: sigma2_icefracturb !--Sigma2 of the ice supersaturation PDF [-]288 REAL, INTENT(OUT), DIMENSION(klon) :: mean_icefracturb !--Mean of the ice supersaturation PDF [-]289 290 REAL, DIMENSION(klon) :: qzero, qsatl, dqsatl, qsati, dqsati !--specific humidity saturation values291 INTEGER :: i292 293 REAL :: qvap_incl, qice_incl, qliq_incl, qiceini_incl !--In-cloud specific quantities [kg/kg]294 REAL :: water_vapor_diff !--Water-vapour diffusion coeff in air (f(T,P)) [m2/s]295 REAL :: air_thermal_conduct !--Thermal conductivity of air (f(T)) [J/m/K/s]296 REAL :: C0 !--Lagrangian structure function [-]297 REAL :: tau_dissipturb298 REAL :: invtau_phaserelax299 REAL :: sigma2_pdf300 REAL :: ai, bi, B0301 REAL :: sursat_iceliq302 REAL :: sursat_equ303 REAL :: liqfra_max304 REAL :: sursat_iceext305 REAL :: nb_crystals !--number concentration of ice crystals [#/m3]306 REAL :: moment1_PSD !--1st moment of ice PSD307 REAL :: N0_PSD, lambda_PSD !--parameters of the exponential PSD308 309 REAL :: cldfra1D310 REAL :: rho_air311 REAL :: psati !--saturation vapor pressure wrt ice [Pa]312 REAL :: sigmaw2 !--variance of vertical turbulent velocity [m2/s2]313 314 REAL :: tempvig1, tempvig2315 316 tempvig1 = -21.06 + RTT317 tempvig2 = -30.35 + RTT318 C0 = 10. !--value assumed in Field2014319 sursat_iceext = -0.1320 qzero(:) = 0.321 cldfraliq(:) = 0.322 qliq(:) = 0.323 qice(:) = 0.324 qvap_cld(:) = 0.325 sigma2_icefracturb(:) = 0.326 mean_icefracturb(:) = 0.327 328 !--wrt liquid329 CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,1,.false.,qsatl(:),dqsatl(:))330 !--wrt ice331 CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,2,.false.,qsati(:),dqsati(:))332 333 334 DO i=1,klon335 rho_air = pplay(i) / temp(i) / RD336 ! assuming turbulence isotropy, tke=3/2*sigmaw2337 sigmaw2=2./3*tke(i)338 ! because cldfra is intent in, but can be locally modified due to test339 cldfra1D = cldfra(i)340 ! activate param for concerned grid points and for cloudy conditions341 IF ((pticefracturb(i)) .AND. (cldfra(i) .GT. 0.)) THEN342 IF (cldfra(i) .GE. 1.0) THEN343 cldfra1D = 1.0344 END IF345 346 ! T>0°C, no ice allowed347 IF ( temp(i) .GE. RTT ) THEN348 qvap_cld(i) = qsatl(i) * cldfra1D349 qliq(i) = MAX(0.0,qtot_incl(i)-qsatl(i)) * cldfra1D350 qice(i) = 0.351 cldfraliq(i) = 1.352 icefrac(i) = 0.353 dicefracdT(i) = 0.354 355 ! T<-38°C, no liquid allowed356 ELSE IF ( temp(i) .LE. temp_nowater) THEN357 qvap_cld(i) = qsati(i) * cldfra1D358 qliq(i) = 0.359 qice(i) = MAX(0.0,qtot_incl(i)-qsati(i)) * cldfra1D360 cldfraliq(i) = 0.361 icefrac(i) = 1.362 dicefracdT(i) = 0.363 364 365 !---------------------------------------------------------366 !-- MIXED PHASE TEMPERATURE REGIME367 !---------------------------------------------------------368 !--In the mixed phase regime (-38°C< T <0°C) we distinguish369 !--3 possible subcases.370 !--1. No pre-existing ice371 !--2A. Pre-existing ice and no turbulence372 !--2B. Pre-existing ice and turbulence373 ELSE374 375 ! gamma_snwretro controls the contribution of snowflakes to the negative feedback376 ! note that for reasons related to inetarctions with the condensation iteration in lscp_main377 ! we consider here the mean snowflake concentration in the mesh (not the in-cloud concentration)378 ! when poprecip is active, it will be worth testing considering the incloud fraction, dividing379 ! by znebprecipcld380 ! qiceini_incl = qice_ini(i) / cldfra1D + &381 ! gamma_snwretro * snowcld(i) * RG * dtime / ( paprsdn(i) - paprsup(i) )382 ! assuming constant snowfall velocity383 qiceini_incl = qice_ini(i) / cldfra1D + gamma_snwretro * snowcld(i) / pplay(i) * RD * temp(i) / snow_fallspeed384 385 !--1. No preexisting ice and no mixing with environment: if vertical motion, fully liquid386 !--cloud else fully iced cloud387 IF ( (qiceini_incl .LT. eps) .AND. (invtau_e(i) .LT. eps) ) THEN388 IF ( (wvel(i)+sqrt(sigmaw2) .GT. eps) .OR. (tke(i) .GT. eps) ) THEN389 qvap_cld(i) = qsatl(i) * cldfra1D390 qliq(i) = MAX(0.,qtot_incl(i)-qsatl(i)) * cldfra1D391 qice(i) = 0.392 cldfraliq(i) = 1.393 icefrac(i) = 0.394 dicefracdT(i) = 0.395 ELSE396 qvap_cld(i) = qsati(i) * cldfra1D397 qliq(i) = 0.398 qice(i) = MAX(0.,qtot_incl(i)-qsati(i)) * cldfra1D399 cldfraliq(i) = 0.400 icefrac(i) = 1.401 dicefracdT(i) = 0.402 ENDIF403 404 405 !--2. Pre-existing ice and/or mixing with environment:computation of ice properties for406 !--feedback407 ELSE408 409 sursat_iceliq = qsatl(i)/qsati(i) - 1.410 psati = qsati(i) * pplay(i) / (RD/RV)411 412 !--We assume an exponential ice PSD whose parameters413 !--are computed following Morrison&Gettelman 2008414 !--Ice number density is assumed equals to INP density415 !--which is for naero5>0 a function of temperature (DeMott 2010)416 !--bi and B0 are microphysical function characterizing417 !--vapor/ice interactions418 !--tau_phase_relax is the typical time of vapor deposition419 !--onto ice crystals420 421 !--For naero5<=0 INP density is derived from the empirical fit422 !--from MARCUS campaign from Vignon 2021423 !--/!\ Note that option is very specific and should be use for424 !--the Southern Ocean and the Antarctic425 426 IF (naero5 .LE. 0) THEN427 IF ( temp(i) .GT. tempvig1 ) THEN428 nb_crystals = 1.e3 * 10**(-0.14*(temp(i)-tempvig1) - 2.88)429 ELSE IF ( temp(i) .GT. tempvig2 ) THEN430 nb_crystals = 1.e3 * 10**(-0.31*(temp(i)-tempvig1) - 2.88)431 ELSE432 nb_crystals = 1.e3 * 10**(0.)433 ENDIF434 ELSE435 nb_crystals = 1.e3 * 5.94e-5 * ( RTT - temp(i) )**3.33 * naero5**(0.0264*(RTT-temp(i))+0.0033)436 ENDIF437 lambda_PSD = ( (RPI*rho_ice*nb_crystals) / (rho_air * MAX(qiceini_incl , eps) ) ) ** (1./3.)438 N0_PSD = nb_crystals * lambda_PSD439 moment1_PSD = N0_PSD/lambda_PSD**2440 441 !--Formulae for air thermal conductivity and water vapor diffusivity442 !--comes respectively from Beard and Pruppacher (1971)443 !--and Hall and Pruppacher (1976)444 445 air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184446 water_vapor_diff = 2.11*1e-5 * ( temp(i) / RTT )**1.94 * ( 101325 / pplay(i) )447 448 bi = 1./((qsati(i)+qsatl(i))/2.) + RLSTT**2 / RCPD / RV / temp(i)**2449 B0 = 4. * RPI * capa_crystal * 1. / ( RLSTT**2 / air_thermal_conduct / RV / temp(i)**2 &450 + RV * temp(i) / psati / water_vapor_diff )451 invtau_phaserelax = bi * B0 * moment1_PSD452 453 ai = RG / RD / temp(i) * ( RD * RLSTT / RCPD / RV / temp(i) - 1. )454 sursat_equ = (ai * wvel(i) + sursat_e(i)*invtau_e(i)) / (invtau_phaserelax + invtau_e(i))455 ! as sursaturation is by definition lower than -1 and456 ! because local supersaturation > 1 are never found in the atmosphere457 458 !--2A. No TKE : stationnary binary solution depending on vertical velocity and mixing with env.459 ! If Sequ > Siw liquid cloud, else ice cloud460 IF ( tke_dissip(i) .LE. eps ) THEN461 sigma2_icefracturb(i)= 0.462 mean_icefracturb(i) = sursat_equ463 IF (sursat_equ .GT. sursat_iceliq) THEN464 qvap_cld(i) = qsatl(i) * cldfra1D465 qliq(i) = MAX(0.,qtot_incl(i)-qsatl(i)) * cldfra1D466 qice(i) = 0.467 cldfraliq(i) = 1.468 icefrac(i) = 0.469 dicefracdT(i) = 0.470 ELSE471 qvap_cld(i) = qsati(i) * cldfra1D472 qliq(i) = 0.473 qice(i) = MAX(0.,qtot_incl(i)-qsati(i)) * cldfra1D474 cldfraliq(i) = 0.475 icefrac(i) = 1.476 dicefracdT(i) = 0.477 ENDIF478 479 !--2B. TKE and ice : ice supersaturation PDF480 !--we compute the cloud liquid properties with a Gaussian PDF481 !--of ice supersaturation F(Si) (Field2014, Furtado2016).482 !--Parameters of the PDF are function of turbulence and483 !--microphysics/existing ice.484 ELSE485 486 !--Tau_dissipturb is the time needed for turbulence to decay487 !--due to viscosity488 tau_dissipturb = gamma_taud * 2. * sigmaw2 / tke_dissip(i) / C0489 490 !--------------------- PDF COMPUTATIONS ---------------------491 !--Formulae for sigma2_pdf (variance), mean of PDF in Raillard2025492 !--cloud liquid fraction and in-cloud liquid content are given493 !--by integrating resp. F(Si) and Si*F(Si)494 !--Liquid is limited by the available water vapor trough a495 !--maximal liquid fraction496 !--qice_ini(i) / cldfra1D = qiceincld without precip497 498 liqfra_max = MAX(0., (MIN (1.,( qtot_incl(i) - (qice_ini(i) / cldfra1D) - qsati(i) * (1 + sursat_iceext ) ) / ( qsatl(i) - qsati(i) ) ) ) )499 sigma2_pdf = 1./2. * ( ai**2 ) * sigmaw2 * tau_dissipturb / (invtau_phaserelax + invtau_e(i))500 ! sursat ranges between -1 and 1, so we prevent sigma2 so exceed 1501 cldfraliq(i) = 0.5 * (1. - erf( ( sursat_iceliq - sursat_equ) / (SQRT(2.* sigma2_pdf) ) ) )502 IF (cldfraliq(i) .GT. liqfra_max) THEN503 cldfraliq(i) = liqfra_max504 ENDIF505 506 qliq_incl = qsati(i) * SQRT(sigma2_pdf) / SQRT(2.*RPI) * EXP( -1.*(sursat_iceliq - sursat_equ)**2. / (2.*sigma2_pdf) ) &507 - qsati(i) * cldfraliq(i) * (sursat_iceliq - sursat_equ )508 509 sigma2_icefracturb(i)= sigma2_pdf510 mean_icefracturb(i) = sursat_equ511 512 !------------ SPECIFIC VAPOR CONTENT AND WATER CONSERVATION ------------513 514 IF ( (qliq_incl .LE. eps) .OR. (cldfraliq(i) .LE. eps) ) THEN515 qliq_incl = 0.516 cldfraliq(i) = 0.517 END IF518 519 !--Specific humidity is the max between qsati and the weighted mean between520 !--qv in MPC patches and qv in ice-only parts. We assume that MPC parts are521 !--always at qsatl and ice-only parts slightly subsaturated (qsati*sursat_iceext+1)522 !--The whole cloud can therefore be supersaturated but never subsaturated.523 524 qvap_incl = MAX(qsati(i), ( 1. - cldfraliq(i) ) * (sursat_iceext + 1.) * qsati(i) + cldfraliq(i) * qsatl(i) )525 526 IF ( qvap_incl .GE. qtot_incl(i) ) THEN527 qvap_incl = qsati(i)528 qliq_incl = MAX(0.0,qtot_incl(i) - qvap_incl)529 qice_incl = 0.530 531 ELSEIF ( (qvap_incl + qliq_incl) .GE. qtot_incl(i) ) THEN532 qliq_incl = MAX(0.0,qtot_incl(i) - qvap_incl)533 qice_incl = 0.534 ELSE535 qice_incl = qtot_incl(i) - qvap_incl - qliq_incl536 END IF537 538 qvap_cld(i) = qvap_incl * cldfra1D539 qliq(i) = qliq_incl * cldfra1D540 qice(i) = qice_incl * cldfra1D541 IF ((qice(i)+qliq(i)) .GT. 0.) THEN542 icefrac(i) = qice(i) / ( qice(i) + qliq(i) )543 ELSE544 icefrac(i) = 1. ! to keep computation of qsat wrt ice in condensation loop in lmdz_lscp_main545 ENDIF546 dicefracdT(i) = 0.547 548 END IF ! Enough TKE549 550 END IF ! End qini551 552 END IF ! ! MPC temperature553 554 END IF ! pticefracturb and cldfra555 556 ENDDO ! klon557 END SUBROUTINE ICEFRAC_LSCP_TURB558 !559 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++560 561 102 562 103 SUBROUTINE CALC_QSAT_ECMWF(klon,temp,qtot,pressure,tref,phase,flagth,qs,dqs) -
LMDZ6/trunk/libf/phylmd/nuage.f90
r5828 r5996 13 13 USE yomcst_mod_h 14 14 USE dimphy 15 USE lmdz_lscp_ tools, only: icefrac_lscp15 USE lmdz_lscp_phase, only: icefrac_lscp 16 16 USE icefrac_lsc_mod ! computes ice fraction (JBM 3/14) 17 17 USE lmdz_lscp_ini, only : iflag_t_glace
Note: See TracChangeset
for help on using the changeset viewer.
