!$gpum horizontal klon 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) END IF 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) END IF 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) END IF IF ((icefrac(i) .EQ. 0) .OR. (icefrac(i) .EQ. 1)) THEN dicefracdT(i) = 0. END IF END IF ! 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. END IF END IF ! 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. END IF END IF ! 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) END IF END IF ! 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. END IF 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) END IF END IF ! 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. END IF END IF END DO ! 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 USE lmdz_lscp_tools, ONLY: calc_qsat_ecmwf 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. END IF !--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.) END IF ELSE nb_crystals = 1.e3*5.94e-5*(RTT - temp(i))**3.33*naero5**(0.0264*(RTT - temp(i)) + 0.0033) END IF 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. END IF !--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 END IF 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 END IF dicefracdT(i) = 0. END IF ! Enough TKE END IF ! End qini END IF ! ! MPC temperature END IF ! pticefracturb and cldfra END DO ! klon END SUBROUTINE ICEFRAC_LSCP_TURB !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ END MODULE lmdz_lscp_phase