Index: LMDZ6/trunk/libf/phylmd/lmdz_call_cloud_optics_prop.f90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lmdz_call_cloud_optics_prop.f90 (revision 5995) +++ LMDZ6/trunk/libf/phylmd/lmdz_call_cloud_optics_prop.f90 (revision 5996) @@ -39,5 +39,5 @@ USE icefrac_lsc_mod ! computes ice fraction (JBM 3/14) - USE lmdz_lscp_tools, only: icefrac_lscp + USE lmdz_lscp_phase, only: icefrac_lscp USE nuage_mod, ONLY: nuage Index: LMDZ6/trunk/libf/phylmd/lmdz_lscp_main.f90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lmdz_lscp_main.f90 (revision 5995) +++ LMDZ6/trunk/libf/phylmd/lmdz_lscp_main.f90 (revision 5996) @@ -71,6 +71,6 @@ ! USE de modules contenant des fonctions. -USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, calc_gammasat -USE lmdz_lscp_tools, ONLY : icefrac_lscp, icefrac_lscp_turb, distance_to_cloud_top +USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, calc_gammasat, distance_to_cloud_top +USE lmdz_lscp_phase, ONLY : icefrac_lscp, icefrac_lscp_turb USE lmdz_lscp_condensation, ONLY : condensation_lognormal, condensation_ice_supersat USE lmdz_lscp_condensation, ONLY : cloudth, cloudth_v3, cloudth_v6, condensation_cloudth Index: LMDZ6/trunk/libf/phylmd/lmdz_lscp_phase.f90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lmdz_lscp_phase.f90 (revision 5996) +++ LMDZ6/trunk/libf/phylmd/lmdz_lscp_phase.f90 (revision 5996) @@ -0,0 +1,474 @@ +MODULE lmdz_lscp_phase + +IMPLICIT NONE +!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +! This module contains routines that determine +! the phase of the clouds computed in lscp +! routines +!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +CONTAINS + +!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + SUBROUTINE ICEFRAC_LSCP(klon, temp, iflag_ice_thermo, distcltop, temp_cltop, icefrac, dicefracdT) +!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + ! Compute the ice fraction (see e.g. Doutriaux-Boucher & Quaas 2004, section 2.2.) + ! as a function of temperature + ! see also Fig 3 of Madeleine et al. 2020, JAMES, 10.1029/2020MS002046 + ! An option is also available to make the cloud phase depend on distance from cloud + ! top. See Lea Raillard's PhD manuscript +!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + + USE print_control_mod, ONLY: lunout, prt_level + USE lmdz_lscp_ini, ONLY: t_glace_min, t_glace_max, exposant_glace, iflag_t_glace + USE lmdz_lscp_ini, ONLY : RTT, dist_liq, temp_nowater + + IMPLICIT NONE + + + INTEGER, INTENT(IN) :: klon !-- number of horizontal grid points + INTEGER, INTENT(IN) :: iflag_ice_thermo !-- option for cloud phase determination + REAL, INTENT(IN), DIMENSION(klon) :: temp !-- temperature [K] + REAL, INTENT(IN), DIMENSION(klon) :: distcltop !-- distance to cloud top [m] + REAL, INTENT(IN), DIMENSION(klon) :: temp_cltop !-- temperature of cloud top [K] + REAL, INTENT(OUT), DIMENSION(klon) :: icefrac !-- ice fraction in clouds [0-1] + REAL, INTENT(OUT), DIMENSION(klon) :: dicefracdT !-- dicefraction/dT [K-1] + + + INTEGER i + REAL liqfrac_tmp, dicefrac_tmp + REAL Dv, denomdep,beta,qsi,dqsidt + LOGICAL ice_thermo + + CHARACTER (len = 20) :: modname = 'lscp_tools' + CHARACTER (len = 80) :: abort_message + + IF ((iflag_t_glace.LT.2)) THEN !.OR. (iflag_t_glace.GT.6)) THEN + abort_message = 'lscp cannot be used if iflag_t_glace<2 or >6' + CALL abort_physic(modname,abort_message,1) + ENDIF + + IF (.NOT.((iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3))) THEN + abort_message = 'lscp cannot be used without ice thermodynamics' + CALL abort_physic(modname,abort_message,1) + ENDIF + + + DO i=1,klon + + ! old function with sole dependence upon temperature + IF (iflag_t_glace .EQ. 2) THEN + liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min) + liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0) + icefrac(i) = (1.0-liqfrac_tmp)**exposant_glace + IF (icefrac(i) .GT.0.) THEN + dicefracdT(i)= exposant_glace * (icefrac(i)**(exposant_glace-1.)) & + / (t_glace_min - t_glace_max) + ENDIF + + IF ((icefrac(i).EQ.0).OR.(icefrac(i).EQ.1)) THEN + dicefracdT(i)=0. + ENDIF + + ENDIF + + ! function of temperature used in CMIP6 physics + IF (iflag_t_glace .EQ. 3) THEN + liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min) + liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0) + icefrac(i) = 1.0-liqfrac_tmp**exposant_glace + IF ((icefrac(i) .GT.0.) .AND. (liqfrac_tmp .GT. 0.)) THEN + dicefracdT(i)= exposant_glace * ((liqfrac_tmp)**(exposant_glace-1.)) & + / (t_glace_min - t_glace_max) + ELSE + dicefracdT(i)=0. + ENDIF + ENDIF + + ! for iflag_t_glace .GE. 4, the liquid fraction depends upon temperature at cloud top + ! and then decreases with decreasing height + + !with linear function of temperature at cloud top + IF (iflag_t_glace .EQ. 4) THEN + liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min) + liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0) + icefrac(i) = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.) + dicefrac_tmp = - temp(i)/(t_glace_max-t_glace_min) + dicefracdT(i) = dicefrac_tmp*exp(-distcltop(i)/dist_liq) + IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN + dicefracdT(i) = 0. + ENDIF + ENDIF + + ! with CMIP6 function of temperature at cloud top + IF ((iflag_t_glace .EQ. 5) .OR. (iflag_t_glace .EQ. 7)) THEN + liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min) + liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0) + liqfrac_tmp = liqfrac_tmp**exposant_glace + icefrac(i) = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.) + IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN + dicefracdT(i) = 0. + ELSE + dicefracdT(i) = exposant_glace*((liqfrac_tmp)**(exposant_glace-1.))/(t_glace_min- t_glace_max) & + *exp(-distcltop(i)/dist_liq) + ENDIF + ENDIF + + ! with modified function of temperature at cloud top + ! to get largere values around 260 K, works well with t_glace_min = 241K + IF (iflag_t_glace .EQ. 6) THEN + IF (temp(i) .GT. t_glace_max) THEN + liqfrac_tmp = 1. + ELSE + liqfrac_tmp = -((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))**2+1. + ENDIF + liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0) + icefrac(i) = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.) + IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN + dicefracdT(i) = 0. + ELSE + dicefracdT(i) = 2*((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))/(t_glace_max-t_glace_min) & + *exp(-distcltop(i)/dist_liq) + ENDIF + ENDIF + + ! if temperature or temperature of cloud top <-40°C, + IF (iflag_t_glace .GE. 4) THEN + IF ((temp_cltop(i) .LE. temp_nowater) .AND. (temp(i) .LE. t_glace_max)) THEN + icefrac(i) = 1. + dicefracdT(i) = 0. + ENDIF + ENDIF + + + ENDDO ! klon + + RETURN + +END SUBROUTINE ICEFRAC_LSCP +!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + +SUBROUTINE ICEFRAC_LSCP_TURB(klon, dtime, pticefracturb, temp, pplay, paprsdn, paprsup, wvel, qice_ini, snowcld, qtot_incl, cldfra, tke, & + tke_dissip, sursat_e, invtau_e, qliq, qvap_cld, qice, icefrac, dicefracdT, cldfraliq, sigma2_icefracturb, mean_icefracturb) +!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + ! Compute the liquid, ice and vapour content (+ice fraction) based + ! on turbulence (see Fields 2014, Furtado 2016, Raillard 2025) + ! L.Raillard (23/09/24) + ! E.Vignon (03/2025) : additional elements for treatment of convective + ! boundary layer clouds + ! References: + ! - Raillard et al. 2026, JAMES, doi:10.22541/essoar.175096287.71557703/v1 + ! - Vignon et al. 2026, ACP, doi:0.5194/egusphere-2025-4641 +!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + USE lmdz_lscp_ini, ONLY : prt_level, lunout + USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI + USE lmdz_lscp_ini, ONLY : seuil_neb, temp_nowater + USE lmdz_lscp_ini, ONLY : naero5, gamma_snwretro, gamma_taud, capa_crystal, rho_ice + USE lmdz_lscp_ini, ONLY : eps, snow_fallspeed + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: klon !--number of horizontal grid points + REAL, INTENT(IN) :: dtime !--time step [s] + LOGICAL, INTENT(IN), DIMENSION(klon) :: pticefracturb !--grid points concerned by this routine + REAL, INTENT(IN), DIMENSION(klon) :: temp !--temperature + REAL, INTENT(IN), DIMENSION(klon) :: pplay !--pressure in the middle of the layer [Pa] + REAL, INTENT(IN), DIMENSION(klon) :: paprsdn !--pressure at the bottom interface of the layer [Pa] + REAL, INTENT(IN), DIMENSION(klon) :: paprsup !--pressure at the top interface of the layer [Pa] + REAL, INTENT(IN), DIMENSION(klon) :: wvel !--vertical velocity [m/s] + REAL, INTENT(IN), DIMENSION(klon) :: qtot_incl !--specific total cloud water in-cloud content [kg/kg] + REAL, INTENT(IN), DIMENSION(klon) :: cldfra !--cloud fraction in gridbox [-] + REAL, INTENT(IN), DIMENSION(klon) :: tke !--turbulent kinetic energy [m2/s2] + REAL, INTENT(IN), DIMENSION(klon) :: tke_dissip !--TKE dissipation [m2/s3] + + REAL, INTENT(IN), DIMENSION(klon) :: qice_ini !--initial specific ice content gridbox-mean [kg/kg] + REAL, INTENT(IN), DIMENSION(klon) :: snowcld !--in-cloud snowfall flux [kg/m2/s] + REAL, INTENT(IN), DIMENSION(klon) :: sursat_e !--environment supersaturation [-] + REAL, INTENT(IN), DIMENSION(klon) :: invtau_e !--inverse time-scale of mixing with environment [s-1] + REAL, INTENT(OUT), DIMENSION(klon) :: qliq !--specific liquid content gridbox-mean [kg/kg] + REAL, INTENT(OUT), DIMENSION(klon) :: qvap_cld !--specific cloud vapor content, gridbox-mean [kg/kg] + REAL, INTENT(OUT), DIMENSION(klon) :: qice !--specific ice content gridbox-mean [kg/kg] + + REAL, INTENT(INOUT), DIMENSION(klon) :: icefrac !--fraction of ice in condensed water [-] + REAL, INTENT(INOUT), DIMENSION(klon) :: dicefracdT + + REAL, INTENT(OUT), DIMENSION(klon) :: cldfraliq !--fraction of cldfra where liquid [-] + REAL, INTENT(OUT), DIMENSION(klon) :: sigma2_icefracturb !--Sigma2 of the ice supersaturation PDF [-] + REAL, INTENT(OUT), DIMENSION(klon) :: mean_icefracturb !--Mean of the ice supersaturation PDF [-] + + REAL, DIMENSION(klon) :: qzero, qsatl, dqsatl, qsati, dqsati !--specific humidity saturation values + INTEGER :: i + + REAL :: qvap_incl, qice_incl, qliq_incl, qiceini_incl !--In-cloud specific quantities [kg/kg] + REAL :: water_vapor_diff !--Water-vapour diffusion coeff in air (f(T,P)) [m2/s] + REAL :: air_thermal_conduct !--Thermal conductivity of air (f(T)) [J/m/K/s] + REAL :: C0 !--Lagrangian structure function [-] + REAL :: tau_dissipturb + REAL :: invtau_phaserelax + REAL :: sigma2_pdf + REAL :: ai, bi, B0 + REAL :: sursat_iceliq + REAL :: sursat_equ + REAL :: liqfra_max + REAL :: sursat_iceext + REAL :: nb_crystals !--number concentration of ice crystals [#/m3] + REAL :: moment1_PSD !--1st moment of ice PSD + REAL :: N0_PSD, lambda_PSD !--parameters of the exponential PSD + + REAL :: cldfra1D + REAL :: rho_air + REAL :: psati !--saturation vapor pressure wrt ice [Pa] + REAL :: sigmaw2 !--variance of vertical turbulent velocity [m2/s2] + + REAL :: tempvig1, tempvig2 + + tempvig1 = -21.06 + RTT + tempvig2 = -30.35 + RTT + C0 = 10. !--value assumed in Field2014 + sursat_iceext = -0.1 + qzero(:) = 0. + cldfraliq(:) = 0. + qliq(:) = 0. + qice(:) = 0. + qvap_cld(:) = 0. + sigma2_icefracturb(:) = 0. + mean_icefracturb(:) = 0. + + !--wrt liquid + CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,1,.false.,qsatl(:),dqsatl(:)) + !--wrt ice + CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,2,.false.,qsati(:),dqsati(:)) + + + DO i=1,klon + rho_air = pplay(i) / temp(i) / RD + ! assuming turbulence isotropy, tke=3/2*sigmaw2 + sigmaw2=2./3*tke(i) + ! because cldfra is intent in, but can be locally modified due to test + cldfra1D = cldfra(i) + ! activate param for concerned grid points and for cloudy conditions + IF ((pticefracturb(i)) .AND. (cldfra(i) .GT. 0.)) THEN + IF (cldfra(i) .GE. 1.0) THEN + cldfra1D = 1.0 + END IF + + ! T>0°C, no ice allowed + IF ( temp(i) .GE. RTT ) THEN + qvap_cld(i) = qsatl(i) * cldfra1D + qliq(i) = MAX(0.0,qtot_incl(i)-qsatl(i)) * cldfra1D + qice(i) = 0. + cldfraliq(i) = 1. + icefrac(i) = 0. + dicefracdT(i) = 0. + + ! T<-38°C, no liquid allowed + ELSE IF ( temp(i) .LE. temp_nowater) THEN + qvap_cld(i) = qsati(i) * cldfra1D + qliq(i) = 0. + qice(i) = MAX(0.0,qtot_incl(i)-qsati(i)) * cldfra1D + cldfraliq(i) = 0. + icefrac(i) = 1. + dicefracdT(i) = 0. + + + !--------------------------------------------------------- + !-- MIXED PHASE TEMPERATURE REGIME + !--------------------------------------------------------- + !--In the mixed phase regime (-38°C< T <0°C) we distinguish + !--3 possible subcases. + !--1. No pre-existing ice + !--2A. Pre-existing ice and no turbulence + !--2B. Pre-existing ice and turbulence + ELSE + + ! gamma_snwretro controls the contribution of snowflakes to the negative feedback + ! note that for reasons related to inetarctions with the condensation iteration in lscp_main + ! we consider here the mean snowflake concentration in the mesh (not the in-cloud concentration) + ! when poprecip is active, it will be worth testing considering the incloud fraction, dividing + ! by znebprecipcld + ! qiceini_incl = qice_ini(i) / cldfra1D + & + ! gamma_snwretro * snowcld(i) * RG * dtime / ( paprsdn(i) - paprsup(i) ) + ! assuming constant snowfall velocity + qiceini_incl = qice_ini(i) / cldfra1D + gamma_snwretro * snowcld(i) / pplay(i) * RD * temp(i) / snow_fallspeed + + !--1. No preexisting ice and no mixing with environment: if vertical motion, fully liquid + !--cloud else fully iced cloud + IF ( (qiceini_incl .LT. eps) .AND. (invtau_e(i) .LT. eps) ) THEN + IF ( (wvel(i)+sqrt(sigmaw2) .GT. eps) .OR. (tke(i) .GT. eps) ) THEN + qvap_cld(i) = qsatl(i) * cldfra1D + qliq(i) = MAX(0.,qtot_incl(i)-qsatl(i)) * cldfra1D + qice(i) = 0. + cldfraliq(i) = 1. + icefrac(i) = 0. + dicefracdT(i) = 0. + ELSE + qvap_cld(i) = qsati(i) * cldfra1D + qliq(i) = 0. + qice(i) = MAX(0.,qtot_incl(i)-qsati(i)) * cldfra1D + cldfraliq(i) = 0. + icefrac(i) = 1. + dicefracdT(i) = 0. + ENDIF + + + !--2. Pre-existing ice and/or mixing with environment:computation of ice properties for + !--feedback + ELSE + + sursat_iceliq = qsatl(i)/qsati(i) - 1. + psati = qsati(i) * pplay(i) / (RD/RV) + + !--We assume an exponential ice PSD whose parameters + !--are computed following Morrison&Gettelman 2008 + !--Ice number density is assumed equals to INP density + !--which is for naero5>0 a function of temperature (DeMott 2010) + !--bi and B0 are microphysical function characterizing + !--vapor/ice interactions + !--tau_phase_relax is the typical time of vapor deposition + !--onto ice crystals + + !--For naero5<=0 INP density is derived from the empirical fit + !--from MARCUS campaign from Vignon 2021 + !--/!\ Note that option is very specific and should be use for + !--the Southern Ocean and the Antarctic + + IF (naero5 .LE. 0) THEN + IF ( temp(i) .GT. tempvig1 ) THEN + nb_crystals = 1.e3 * 10**(-0.14*(temp(i)-tempvig1) - 2.88) + ELSE IF ( temp(i) .GT. tempvig2 ) THEN + nb_crystals = 1.e3 * 10**(-0.31*(temp(i)-tempvig1) - 2.88) + ELSE + nb_crystals = 1.e3 * 10**(0.) + ENDIF + ELSE + nb_crystals = 1.e3 * 5.94e-5 * ( RTT - temp(i) )**3.33 * naero5**(0.0264*(RTT-temp(i))+0.0033) + ENDIF + lambda_PSD = ( (RPI*rho_ice*nb_crystals) / (rho_air * MAX(qiceini_incl , eps) ) ) ** (1./3.) + N0_PSD = nb_crystals * lambda_PSD + moment1_PSD = N0_PSD/lambda_PSD**2 + + !--Formulae for air thermal conductivity and water vapor diffusivity + !--comes respectively from Beard and Pruppacher (1971) + !--and Hall and Pruppacher (1976) + + air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184 + water_vapor_diff = 2.11*1e-5 * ( temp(i) / RTT )**1.94 * ( 101325 / pplay(i) ) + + bi = 1./((qsati(i)+qsatl(i))/2.) + RLSTT**2 / RCPD / RV / temp(i)**2 + B0 = 4. * RPI * capa_crystal * 1. / ( RLSTT**2 / air_thermal_conduct / RV / temp(i)**2 & + + RV * temp(i) / psati / water_vapor_diff ) + invtau_phaserelax = bi * B0 * moment1_PSD + + ai = RG / RD / temp(i) * ( RD * RLSTT / RCPD / RV / temp(i) - 1. ) + sursat_equ = (ai * wvel(i) + sursat_e(i)*invtau_e(i)) / (invtau_phaserelax + invtau_e(i)) + ! as sursaturation is by definition lower than -1 and + ! because local supersaturation > 1 are never found in the atmosphere + + !--2A. No TKE : stationnary binary solution depending on vertical velocity and mixing with env. + ! If Sequ > Siw liquid cloud, else ice cloud + IF ( tke_dissip(i) .LE. eps ) THEN + sigma2_icefracturb(i)= 0. + mean_icefracturb(i) = sursat_equ + IF (sursat_equ .GT. sursat_iceliq) THEN + qvap_cld(i) = qsatl(i) * cldfra1D + qliq(i) = MAX(0.,qtot_incl(i)-qsatl(i)) * cldfra1D + qice(i) = 0. + cldfraliq(i) = 1. + icefrac(i) = 0. + dicefracdT(i) = 0. + ELSE + qvap_cld(i) = qsati(i) * cldfra1D + qliq(i) = 0. + qice(i) = MAX(0.,qtot_incl(i)-qsati(i)) * cldfra1D + cldfraliq(i) = 0. + icefrac(i) = 1. + dicefracdT(i) = 0. + ENDIF + + !--2B. TKE and ice : ice supersaturation PDF + !--we compute the cloud liquid properties with a Gaussian PDF + !--of ice supersaturation F(Si) (Field2014, Furtado2016). + !--Parameters of the PDF are function of turbulence and + !--microphysics/existing ice. + ELSE + + !--Tau_dissipturb is the time needed for turbulence to decay + !--due to viscosity + tau_dissipturb = gamma_taud * 2. * sigmaw2 / tke_dissip(i) / C0 + + !--------------------- PDF COMPUTATIONS --------------------- + !--Formulae for sigma2_pdf (variance), mean of PDF in Raillard2025 + !--cloud liquid fraction and in-cloud liquid content are given + !--by integrating resp. F(Si) and Si*F(Si) + !--Liquid is limited by the available water vapor trough a + !--maximal liquid fraction + !--qice_ini(i) / cldfra1D = qiceincld without precip + + liqfra_max = MAX(0., (MIN (1.,( qtot_incl(i) - (qice_ini(i) / cldfra1D) - qsati(i) * (1 + sursat_iceext ) ) / ( qsatl(i) - qsati(i) ) ) ) ) + sigma2_pdf = 1./2. * ( ai**2 ) * sigmaw2 * tau_dissipturb / (invtau_phaserelax + invtau_e(i)) + ! sursat ranges between -1 and 1, so we prevent sigma2 so exceed 1 + cldfraliq(i) = 0.5 * (1. - erf( ( sursat_iceliq - sursat_equ) / (SQRT(2.* sigma2_pdf) ) ) ) + IF (cldfraliq(i) .GT. liqfra_max) THEN + cldfraliq(i) = liqfra_max + ENDIF + + qliq_incl = qsati(i) * SQRT(sigma2_pdf) / SQRT(2.*RPI) * EXP( -1.*(sursat_iceliq - sursat_equ)**2. / (2.*sigma2_pdf) ) & + - qsati(i) * cldfraliq(i) * (sursat_iceliq - sursat_equ ) + + sigma2_icefracturb(i)= sigma2_pdf + mean_icefracturb(i) = sursat_equ + + !------------ SPECIFIC VAPOR CONTENT AND WATER CONSERVATION ------------ + + IF ( (qliq_incl .LE. eps) .OR. (cldfraliq(i) .LE. eps) ) THEN + qliq_incl = 0. + cldfraliq(i) = 0. + END IF + + !--Specific humidity is the max between qsati and the weighted mean between + !--qv in MPC patches and qv in ice-only parts. We assume that MPC parts are + !--always at qsatl and ice-only parts slightly subsaturated (qsati*sursat_iceext+1) + !--The whole cloud can therefore be supersaturated but never subsaturated. + + qvap_incl = MAX(qsati(i), ( 1. - cldfraliq(i) ) * (sursat_iceext + 1.) * qsati(i) + cldfraliq(i) * qsatl(i) ) + + IF ( qvap_incl .GE. qtot_incl(i) ) THEN + qvap_incl = qsati(i) + qliq_incl = MAX(0.0,qtot_incl(i) - qvap_incl) + qice_incl = 0. + + ELSEIF ( (qvap_incl + qliq_incl) .GE. qtot_incl(i) ) THEN + qliq_incl = MAX(0.0,qtot_incl(i) - qvap_incl) + qice_incl = 0. + ELSE + qice_incl = qtot_incl(i) - qvap_incl - qliq_incl + END IF + + qvap_cld(i) = qvap_incl * cldfra1D + qliq(i) = qliq_incl * cldfra1D + qice(i) = qice_incl * cldfra1D + IF ((qice(i)+qliq(i)) .GT. 0.) THEN + icefrac(i) = qice(i) / ( qice(i) + qliq(i) ) + ELSE + icefrac(i) = 1. ! to keep computation of qsat wrt ice in condensation loop in lmdz_lscp_main + ENDIF + dicefracdT(i) = 0. + + END IF ! Enough TKE + + END IF ! End qini + + END IF ! ! MPC temperature + + END IF ! pticefracturb and cldfra + + ENDDO ! klon + +END SUBROUTINE ICEFRAC_LSCP_TURB +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + +END MODULE lmdz_lscp_phase Index: LMDZ6/trunk/libf/phylmd/lmdz_lscp_tools.f90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lmdz_lscp_tools.f90 (revision 5995) +++ LMDZ6/trunk/libf/phylmd/lmdz_lscp_tools.f90 (revision 5996) @@ -100,463 +100,4 @@ END SUBROUTINE FALLICE_VELOCITY !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -SUBROUTINE ICEFRAC_LSCP(klon, temp, iflag_ice_thermo, distcltop, temp_cltop, icefrac, dicefracdT) -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - ! Compute the ice fraction 1-xliq (see e.g. - ! Doutriaux-Boucher & Quaas 2004, section 2.2.) - ! as a function of temperature - ! see also Fig 3 of Madeleine et al. 2020, JAMES -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - - USE print_control_mod, ONLY: lunout, prt_level - USE lmdz_lscp_ini, ONLY: t_glace_min, t_glace_max, exposant_glace, iflag_t_glace - USE lmdz_lscp_ini, ONLY : RTT, dist_liq, temp_nowater - - IMPLICIT NONE - - - INTEGER, INTENT(IN) :: klon ! number of horizontal grid points - REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature - REAL, INTENT(IN), DIMENSION(klon) :: distcltop ! distance to cloud top - REAL, INTENT(IN), DIMENSION(klon) :: temp_cltop ! temperature of cloud top - INTEGER, INTENT(IN) :: iflag_ice_thermo - REAL, INTENT(OUT), DIMENSION(klon) :: icefrac - REAL, INTENT(OUT), DIMENSION(klon) :: dicefracdT - - - INTEGER i - REAL liqfrac_tmp, dicefrac_tmp - REAL Dv, denomdep,beta,qsi,dqsidt - LOGICAL ice_thermo - - CHARACTER (len = 20) :: modname = 'lscp_tools' - CHARACTER (len = 80) :: abort_message - - IF ((iflag_t_glace.LT.2)) THEN !.OR. (iflag_t_glace.GT.6)) THEN - abort_message = 'lscp cannot be used if iflag_t_glace<2 or >6' - CALL abort_physic(modname,abort_message,1) - ENDIF - - IF (.NOT.((iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3))) THEN - abort_message = 'lscp cannot be used without ice thermodynamics' - CALL abort_physic(modname,abort_message,1) - ENDIF - - - DO i=1,klon - - ! old function with sole dependence upon temperature - IF (iflag_t_glace .EQ. 2) THEN - liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min) - liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0) - icefrac(i) = (1.0-liqfrac_tmp)**exposant_glace - IF (icefrac(i) .GT.0.) THEN - dicefracdT(i)= exposant_glace * (icefrac(i)**(exposant_glace-1.)) & - / (t_glace_min - t_glace_max) - ENDIF - - IF ((icefrac(i).EQ.0).OR.(icefrac(i).EQ.1)) THEN - dicefracdT(i)=0. - ENDIF - - ENDIF - - ! function of temperature used in CMIP6 physics - IF (iflag_t_glace .EQ. 3) THEN - liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min) - liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0) - icefrac(i) = 1.0-liqfrac_tmp**exposant_glace - IF ((icefrac(i) .GT.0.) .AND. (liqfrac_tmp .GT. 0.)) THEN - dicefracdT(i)= exposant_glace * ((liqfrac_tmp)**(exposant_glace-1.)) & - / (t_glace_min - t_glace_max) - ELSE - dicefracdT(i)=0. - ENDIF - ENDIF - - ! for iflag_t_glace .GE. 4, the liquid fraction depends upon temperature at cloud top - ! and then decreases with decreasing height - - !with linear function of temperature at cloud top - IF (iflag_t_glace .EQ. 4) THEN - liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min) - liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0) - icefrac(i) = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.) - dicefrac_tmp = - temp(i)/(t_glace_max-t_glace_min) - dicefracdT(i) = dicefrac_tmp*exp(-distcltop(i)/dist_liq) - IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN - dicefracdT(i) = 0. - ENDIF - ENDIF - - ! with CMIP6 function of temperature at cloud top - IF ((iflag_t_glace .EQ. 5) .OR. (iflag_t_glace .EQ. 7)) THEN - liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min) - liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0) - liqfrac_tmp = liqfrac_tmp**exposant_glace - icefrac(i) = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.) - IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN - dicefracdT(i) = 0. - ELSE - dicefracdT(i) = exposant_glace*((liqfrac_tmp)**(exposant_glace-1.))/(t_glace_min- t_glace_max) & - *exp(-distcltop(i)/dist_liq) - ENDIF - ENDIF - - ! with modified function of temperature at cloud top - ! to get largere values around 260 K, works well with t_glace_min = 241K - IF (iflag_t_glace .EQ. 6) THEN - IF (temp(i) .GT. t_glace_max) THEN - liqfrac_tmp = 1. - ELSE - liqfrac_tmp = -((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))**2+1. - ENDIF - liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0) - icefrac(i) = MAX(MIN(1.,1.0 - liqfrac_tmp*exp(-distcltop(i)/dist_liq)),0.) - IF ((liqfrac_tmp .LE.0) .OR. (liqfrac_tmp .GE. 1)) THEN - dicefracdT(i) = 0. - ELSE - dicefracdT(i) = 2*((temp(i)-t_glace_max) / (t_glace_max-t_glace_min))/(t_glace_max-t_glace_min) & - *exp(-distcltop(i)/dist_liq) - ENDIF - ENDIF - - ! if temperature or temperature of cloud top <-40°C, - IF (iflag_t_glace .GE. 4) THEN - IF ((temp_cltop(i) .LE. temp_nowater) .AND. (temp(i) .LE. t_glace_max)) THEN - icefrac(i) = 1. - dicefracdT(i) = 0. - ENDIF - ENDIF - - - ENDDO ! klon - RETURN - -END SUBROUTINE ICEFRAC_LSCP -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - -SUBROUTINE ICEFRAC_LSCP_TURB(klon, dtime, pticefracturb, temp, pplay, paprsdn, paprsup, wvel, qice_ini, snowcld, qtot_incl, cldfra, tke, & - tke_dissip, sursat_e, invtau_e, qliq, qvap_cld, qice, icefrac, dicefracdT, cldfraliq, sigma2_icefracturb, mean_icefracturb) -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - ! Compute the liquid, ice and vapour content (+ice fraction) based - ! on turbulence (see Fields 2014, Furtado 2016, Raillard 2025) - ! L.Raillard (23/09/24) - ! E.Vignon (03/2025) : additional elements for treatment of convective - ! boundary layer clouds -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - - - USE lmdz_lscp_ini, ONLY : prt_level, lunout - USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI - USE lmdz_lscp_ini, ONLY : seuil_neb, temp_nowater - USE lmdz_lscp_ini, ONLY : naero5, gamma_snwretro, gamma_taud, capa_crystal, rho_ice - USE lmdz_lscp_ini, ONLY : eps, snow_fallspeed - - IMPLICIT NONE - - INTEGER, INTENT(IN) :: klon !--number of horizontal grid points - REAL, INTENT(IN) :: dtime !--time step [s] - LOGICAL, INTENT(IN), DIMENSION(klon) :: pticefracturb !--grid points concerned by this routine - REAL, INTENT(IN), DIMENSION(klon) :: temp !--temperature - REAL, INTENT(IN), DIMENSION(klon) :: pplay !--pressure in the middle of the layer [Pa] - REAL, INTENT(IN), DIMENSION(klon) :: paprsdn !--pressure at the bottom interface of the layer [Pa] - REAL, INTENT(IN), DIMENSION(klon) :: paprsup !--pressure at the top interface of the layer [Pa] - REAL, INTENT(IN), DIMENSION(klon) :: wvel !--vertical velocity [m/s] - REAL, INTENT(IN), DIMENSION(klon) :: qtot_incl !--specific total cloud water in-cloud content [kg/kg] - REAL, INTENT(IN), DIMENSION(klon) :: cldfra !--cloud fraction in gridbox [-] - REAL, INTENT(IN), DIMENSION(klon) :: tke !--turbulent kinetic energy [m2/s2] - REAL, INTENT(IN), DIMENSION(klon) :: tke_dissip !--TKE dissipation [m2/s3] - - REAL, INTENT(IN), DIMENSION(klon) :: qice_ini !--initial specific ice content gridbox-mean [kg/kg] - REAL, INTENT(IN), DIMENSION(klon) :: snowcld !--in-cloud snowfall flux [kg/m2/s] - REAL, INTENT(IN), DIMENSION(klon) :: sursat_e !--environment supersaturation [-] - REAL, INTENT(IN), DIMENSION(klon) :: invtau_e !--inverse time-scale of mixing with environment [s-1] - REAL, INTENT(OUT), DIMENSION(klon) :: qliq !--specific liquid content gridbox-mean [kg/kg] - REAL, INTENT(OUT), DIMENSION(klon) :: qvap_cld !--specific cloud vapor content, gridbox-mean [kg/kg] - REAL, INTENT(OUT), DIMENSION(klon) :: qice !--specific ice content gridbox-mean [kg/kg] - - REAL, INTENT(INOUT), DIMENSION(klon) :: icefrac !--fraction of ice in condensed water [-] - REAL, INTENT(INOUT), DIMENSION(klon) :: dicefracdT - - REAL, INTENT(OUT), DIMENSION(klon) :: cldfraliq !--fraction of cldfra where liquid [-] - REAL, INTENT(OUT), DIMENSION(klon) :: sigma2_icefracturb !--Sigma2 of the ice supersaturation PDF [-] - REAL, INTENT(OUT), DIMENSION(klon) :: mean_icefracturb !--Mean of the ice supersaturation PDF [-] - - REAL, DIMENSION(klon) :: qzero, qsatl, dqsatl, qsati, dqsati !--specific humidity saturation values - INTEGER :: i - - REAL :: qvap_incl, qice_incl, qliq_incl, qiceini_incl !--In-cloud specific quantities [kg/kg] - REAL :: water_vapor_diff !--Water-vapour diffusion coeff in air (f(T,P)) [m2/s] - REAL :: air_thermal_conduct !--Thermal conductivity of air (f(T)) [J/m/K/s] - REAL :: C0 !--Lagrangian structure function [-] - REAL :: tau_dissipturb - REAL :: invtau_phaserelax - REAL :: sigma2_pdf - REAL :: ai, bi, B0 - REAL :: sursat_iceliq - REAL :: sursat_equ - REAL :: liqfra_max - REAL :: sursat_iceext - REAL :: nb_crystals !--number concentration of ice crystals [#/m3] - REAL :: moment1_PSD !--1st moment of ice PSD - REAL :: N0_PSD, lambda_PSD !--parameters of the exponential PSD - - REAL :: cldfra1D - REAL :: rho_air - REAL :: psati !--saturation vapor pressure wrt ice [Pa] - REAL :: sigmaw2 !--variance of vertical turbulent velocity [m2/s2] - - REAL :: tempvig1, tempvig2 - - tempvig1 = -21.06 + RTT - tempvig2 = -30.35 + RTT - C0 = 10. !--value assumed in Field2014 - sursat_iceext = -0.1 - qzero(:) = 0. - cldfraliq(:) = 0. - qliq(:) = 0. - qice(:) = 0. - qvap_cld(:) = 0. - sigma2_icefracturb(:) = 0. - mean_icefracturb(:) = 0. - - !--wrt liquid - CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,1,.false.,qsatl(:),dqsatl(:)) - !--wrt ice - CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,2,.false.,qsati(:),dqsati(:)) - - - DO i=1,klon - rho_air = pplay(i) / temp(i) / RD - ! assuming turbulence isotropy, tke=3/2*sigmaw2 - sigmaw2=2./3*tke(i) - ! because cldfra is intent in, but can be locally modified due to test - cldfra1D = cldfra(i) - ! activate param for concerned grid points and for cloudy conditions - IF ((pticefracturb(i)) .AND. (cldfra(i) .GT. 0.)) THEN - IF (cldfra(i) .GE. 1.0) THEN - cldfra1D = 1.0 - END IF - - ! T>0°C, no ice allowed - IF ( temp(i) .GE. RTT ) THEN - qvap_cld(i) = qsatl(i) * cldfra1D - qliq(i) = MAX(0.0,qtot_incl(i)-qsatl(i)) * cldfra1D - qice(i) = 0. - cldfraliq(i) = 1. - icefrac(i) = 0. - dicefracdT(i) = 0. - - ! T<-38°C, no liquid allowed - ELSE IF ( temp(i) .LE. temp_nowater) THEN - qvap_cld(i) = qsati(i) * cldfra1D - qliq(i) = 0. - qice(i) = MAX(0.0,qtot_incl(i)-qsati(i)) * cldfra1D - cldfraliq(i) = 0. - icefrac(i) = 1. - dicefracdT(i) = 0. - - - !--------------------------------------------------------- - !-- MIXED PHASE TEMPERATURE REGIME - !--------------------------------------------------------- - !--In the mixed phase regime (-38°C< T <0°C) we distinguish - !--3 possible subcases. - !--1. No pre-existing ice - !--2A. Pre-existing ice and no turbulence - !--2B. Pre-existing ice and turbulence - ELSE - - ! gamma_snwretro controls the contribution of snowflakes to the negative feedback - ! note that for reasons related to inetarctions with the condensation iteration in lscp_main - ! we consider here the mean snowflake concentration in the mesh (not the in-cloud concentration) - ! when poprecip is active, it will be worth testing considering the incloud fraction, dividing - ! by znebprecipcld - ! qiceini_incl = qice_ini(i) / cldfra1D + & - ! gamma_snwretro * snowcld(i) * RG * dtime / ( paprsdn(i) - paprsup(i) ) - ! assuming constant snowfall velocity - qiceini_incl = qice_ini(i) / cldfra1D + gamma_snwretro * snowcld(i) / pplay(i) * RD * temp(i) / snow_fallspeed - - !--1. No preexisting ice and no mixing with environment: if vertical motion, fully liquid - !--cloud else fully iced cloud - IF ( (qiceini_incl .LT. eps) .AND. (invtau_e(i) .LT. eps) ) THEN - IF ( (wvel(i)+sqrt(sigmaw2) .GT. eps) .OR. (tke(i) .GT. eps) ) THEN - qvap_cld(i) = qsatl(i) * cldfra1D - qliq(i) = MAX(0.,qtot_incl(i)-qsatl(i)) * cldfra1D - qice(i) = 0. - cldfraliq(i) = 1. - icefrac(i) = 0. - dicefracdT(i) = 0. - ELSE - qvap_cld(i) = qsati(i) * cldfra1D - qliq(i) = 0. - qice(i) = MAX(0.,qtot_incl(i)-qsati(i)) * cldfra1D - cldfraliq(i) = 0. - icefrac(i) = 1. - dicefracdT(i) = 0. - ENDIF - - - !--2. Pre-existing ice and/or mixing with environment:computation of ice properties for - !--feedback - ELSE - - sursat_iceliq = qsatl(i)/qsati(i) - 1. - psati = qsati(i) * pplay(i) / (RD/RV) - - !--We assume an exponential ice PSD whose parameters - !--are computed following Morrison&Gettelman 2008 - !--Ice number density is assumed equals to INP density - !--which is for naero5>0 a function of temperature (DeMott 2010) - !--bi and B0 are microphysical function characterizing - !--vapor/ice interactions - !--tau_phase_relax is the typical time of vapor deposition - !--onto ice crystals - - !--For naero5<=0 INP density is derived from the empirical fit - !--from MARCUS campaign from Vignon 2021 - !--/!\ Note that option is very specific and should be use for - !--the Southern Ocean and the Antarctic - - IF (naero5 .LE. 0) THEN - IF ( temp(i) .GT. tempvig1 ) THEN - nb_crystals = 1.e3 * 10**(-0.14*(temp(i)-tempvig1) - 2.88) - ELSE IF ( temp(i) .GT. tempvig2 ) THEN - nb_crystals = 1.e3 * 10**(-0.31*(temp(i)-tempvig1) - 2.88) - ELSE - nb_crystals = 1.e3 * 10**(0.) - ENDIF - ELSE - nb_crystals = 1.e3 * 5.94e-5 * ( RTT - temp(i) )**3.33 * naero5**(0.0264*(RTT-temp(i))+0.0033) - ENDIF - lambda_PSD = ( (RPI*rho_ice*nb_crystals) / (rho_air * MAX(qiceini_incl , eps) ) ) ** (1./3.) - N0_PSD = nb_crystals * lambda_PSD - moment1_PSD = N0_PSD/lambda_PSD**2 - - !--Formulae for air thermal conductivity and water vapor diffusivity - !--comes respectively from Beard and Pruppacher (1971) - !--and Hall and Pruppacher (1976) - - air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184 - water_vapor_diff = 2.11*1e-5 * ( temp(i) / RTT )**1.94 * ( 101325 / pplay(i) ) - - bi = 1./((qsati(i)+qsatl(i))/2.) + RLSTT**2 / RCPD / RV / temp(i)**2 - B0 = 4. * RPI * capa_crystal * 1. / ( RLSTT**2 / air_thermal_conduct / RV / temp(i)**2 & - + RV * temp(i) / psati / water_vapor_diff ) - invtau_phaserelax = bi * B0 * moment1_PSD - - ai = RG / RD / temp(i) * ( RD * RLSTT / RCPD / RV / temp(i) - 1. ) - sursat_equ = (ai * wvel(i) + sursat_e(i)*invtau_e(i)) / (invtau_phaserelax + invtau_e(i)) - ! as sursaturation is by definition lower than -1 and - ! because local supersaturation > 1 are never found in the atmosphere - - !--2A. No TKE : stationnary binary solution depending on vertical velocity and mixing with env. - ! If Sequ > Siw liquid cloud, else ice cloud - IF ( tke_dissip(i) .LE. eps ) THEN - sigma2_icefracturb(i)= 0. - mean_icefracturb(i) = sursat_equ - IF (sursat_equ .GT. sursat_iceliq) THEN - qvap_cld(i) = qsatl(i) * cldfra1D - qliq(i) = MAX(0.,qtot_incl(i)-qsatl(i)) * cldfra1D - qice(i) = 0. - cldfraliq(i) = 1. - icefrac(i) = 0. - dicefracdT(i) = 0. - ELSE - qvap_cld(i) = qsati(i) * cldfra1D - qliq(i) = 0. - qice(i) = MAX(0.,qtot_incl(i)-qsati(i)) * cldfra1D - cldfraliq(i) = 0. - icefrac(i) = 1. - dicefracdT(i) = 0. - ENDIF - - !--2B. TKE and ice : ice supersaturation PDF - !--we compute the cloud liquid properties with a Gaussian PDF - !--of ice supersaturation F(Si) (Field2014, Furtado2016). - !--Parameters of the PDF are function of turbulence and - !--microphysics/existing ice. - ELSE - - !--Tau_dissipturb is the time needed for turbulence to decay - !--due to viscosity - tau_dissipturb = gamma_taud * 2. * sigmaw2 / tke_dissip(i) / C0 - - !--------------------- PDF COMPUTATIONS --------------------- - !--Formulae for sigma2_pdf (variance), mean of PDF in Raillard2025 - !--cloud liquid fraction and in-cloud liquid content are given - !--by integrating resp. F(Si) and Si*F(Si) - !--Liquid is limited by the available water vapor trough a - !--maximal liquid fraction - !--qice_ini(i) / cldfra1D = qiceincld without precip - - liqfra_max = MAX(0., (MIN (1.,( qtot_incl(i) - (qice_ini(i) / cldfra1D) - qsati(i) * (1 + sursat_iceext ) ) / ( qsatl(i) - qsati(i) ) ) ) ) - sigma2_pdf = 1./2. * ( ai**2 ) * sigmaw2 * tau_dissipturb / (invtau_phaserelax + invtau_e(i)) - ! sursat ranges between -1 and 1, so we prevent sigma2 so exceed 1 - cldfraliq(i) = 0.5 * (1. - erf( ( sursat_iceliq - sursat_equ) / (SQRT(2.* sigma2_pdf) ) ) ) - IF (cldfraliq(i) .GT. liqfra_max) THEN - cldfraliq(i) = liqfra_max - ENDIF - - qliq_incl = qsati(i) * SQRT(sigma2_pdf) / SQRT(2.*RPI) * EXP( -1.*(sursat_iceliq - sursat_equ)**2. / (2.*sigma2_pdf) ) & - - qsati(i) * cldfraliq(i) * (sursat_iceliq - sursat_equ ) - - sigma2_icefracturb(i)= sigma2_pdf - mean_icefracturb(i) = sursat_equ - - !------------ SPECIFIC VAPOR CONTENT AND WATER CONSERVATION ------------ - - IF ( (qliq_incl .LE. eps) .OR. (cldfraliq(i) .LE. eps) ) THEN - qliq_incl = 0. - cldfraliq(i) = 0. - END IF - - !--Specific humidity is the max between qsati and the weighted mean between - !--qv in MPC patches and qv in ice-only parts. We assume that MPC parts are - !--always at qsatl and ice-only parts slightly subsaturated (qsati*sursat_iceext+1) - !--The whole cloud can therefore be supersaturated but never subsaturated. - - qvap_incl = MAX(qsati(i), ( 1. - cldfraliq(i) ) * (sursat_iceext + 1.) * qsati(i) + cldfraliq(i) * qsatl(i) ) - - IF ( qvap_incl .GE. qtot_incl(i) ) THEN - qvap_incl = qsati(i) - qliq_incl = MAX(0.0,qtot_incl(i) - qvap_incl) - qice_incl = 0. - - ELSEIF ( (qvap_incl + qliq_incl) .GE. qtot_incl(i) ) THEN - qliq_incl = MAX(0.0,qtot_incl(i) - qvap_incl) - qice_incl = 0. - ELSE - qice_incl = qtot_incl(i) - qvap_incl - qliq_incl - END IF - - qvap_cld(i) = qvap_incl * cldfra1D - qliq(i) = qliq_incl * cldfra1D - qice(i) = qice_incl * cldfra1D - IF ((qice(i)+qliq(i)) .GT. 0.) THEN - icefrac(i) = qice(i) / ( qice(i) + qliq(i) ) - ELSE - icefrac(i) = 1. ! to keep computation of qsat wrt ice in condensation loop in lmdz_lscp_main - ENDIF - dicefracdT(i) = 0. - - END IF ! Enough TKE - - END IF ! End qini - - END IF ! ! MPC temperature - - END IF ! pticefracturb and cldfra - - ENDDO ! klon -END SUBROUTINE ICEFRAC_LSCP_TURB -! -!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - SUBROUTINE CALC_QSAT_ECMWF(klon,temp,qtot,pressure,tref,phase,flagth,qs,dqs) Index: LMDZ6/trunk/libf/phylmd/nuage.f90 =================================================================== --- LMDZ6/trunk/libf/phylmd/nuage.f90 (revision 5995) +++ LMDZ6/trunk/libf/phylmd/nuage.f90 (revision 5996) @@ -13,5 +13,5 @@ USE yomcst_mod_h USE dimphy - USE lmdz_lscp_tools, only: icefrac_lscp + USE lmdz_lscp_phase, only: icefrac_lscp USE icefrac_lsc_mod ! computes ice fraction (JBM 3/14) USE lmdz_lscp_ini, only : iflag_t_glace Index: LMDZ6/trunk/libf/phylmdiso/lmdz_lscp_phase.f90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/lmdz_lscp_phase.f90 (revision 5996) +++ LMDZ6/trunk/libf/phylmdiso/lmdz_lscp_phase.f90 (revision 5996) @@ -0,0 +1,1 @@ +link ../phylmd/lmdz_lscp_phase.f90