MODULE lmdz_lscp_tools IMPLICIT NONE CONTAINS !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE FALLICE_VELOCITY(klon,iwc,temp,rho,pres,ptconv,velo) ! Ref: ! Stubenrauch, C. J., Bonazzola, M., ! Protopapadaki, S. E., & Musat, I. (2019). ! New cloud system metrics to assess bulk ! ice cloud schemes in a GCM. Journal of ! Advances in Modeling Earth Systems, 11, ! 3212–3234. https://doi.org/10.1029/2019MS001642 use lmdz_lscp_ini, only: iflag_vice, ffallv_con, ffallv_lsc use lmdz_lscp_ini, only: cice_velo, dice_velo IMPLICIT NONE INTEGER, INTENT(IN) :: klon REAL, INTENT(IN), DIMENSION(klon) :: iwc ! specific ice water content [kg/m3] REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature [K] REAL, INTENT(IN), DIMENSION(klon) :: rho ! dry air density [kg/m3] REAL, INTENT(IN), DIMENSION(klon) :: pres ! air pressure [Pa] LOGICAL, INTENT(IN), DIMENSION(klon) :: ptconv ! convective point [-] REAL, INTENT(OUT), DIMENSION(klon) :: velo ! fallspeed velocity of crystals [m/s] INTEGER i REAL logvm,iwcg,tempc,phpa,fallv_tun REAL m2ice, m2snow, vmice, vmsnow REAL aice, bice, asnow, bsnow DO i=1,klon IF (ptconv(i)) THEN fallv_tun=ffallv_con ELSE fallv_tun=ffallv_lsc ENDIF tempc=temp(i)-273.15 ! celcius temp iwcg=MAX(iwc(i)*1000.,1E-3) ! iwc in g/m3. We set a minimum value to prevent from division by 0 phpa=pres(i)/100. ! pressure in hPa IF (iflag_vice .EQ. 1) THEN ! so-called 'empirical parameterization' in Stubenrauch et al. 2019 if (tempc .GE. -60.0) then logvm= -0.0000414122*tempc*tempc*log(iwcg)-0.00538922*tempc*log(iwcg) & -0.0516344*log(iwcg)+0.00216078*tempc + 1.9714 velo(i)=exp(logvm) else velo(i)=65.0*(iwcg**0.2)*(150./phpa)**0.15 endif velo(i)=fallv_tun*velo(i)/100.0 ! from cm/s to m/s ELSE IF (iflag_vice .EQ. 2) THEN ! so called PSDM empirical coherent bulk ice scheme in Stubenrauch et al. 2019 aice=0.587 bice=2.45 asnow=0.0444 bsnow=2.1 m2ice=((iwcg*0.001/aice)/(exp(13.6-bice*7.76+0.479*bice**2)* & exp((-0.0361+bice*0.0151+0.00149*bice**2)*tempc))) & **(1./(0.807+bice*0.00581+0.0457*bice**2)) vmice=100.*1042.4*exp(13.6-(bice+1)*7.76+0.479*(bice+1.)**2)*exp((-0.0361+ & (bice+1.)*0.0151+0.00149*(bice+1.)**2)*tempc) & *(m2ice**(0.807+(bice+1.)*0.00581+0.0457*(bice+1.)**2))/(iwcg*0.001/aice) vmice=vmice*((1000./phpa)**0.2) m2snow=((iwcg*0.001/asnow)/(exp(13.6-bsnow*7.76+0.479*bsnow**2)* & exp((-0.0361+bsnow*0.0151+0.00149*bsnow**2)*tempc))) & **(1./(0.807+bsnow*0.00581+0.0457*bsnow**2)) vmsnow=100.*14.3*exp(13.6-(bsnow+.416)*7.76+0.479*(bsnow+.416)**2)& *exp((-0.0361+(bsnow+.416)*0.0151+0.00149*(bsnow+.416)**2)*tempc)& *(m2snow**(0.807+(bsnow+.416)*0.00581+0.0457*(bsnow+.416)**2))/(iwcg*0.001/asnow) vmsnow=vmsnow*((1000./phpa)**0.35) velo(i)=fallv_tun*min(vmsnow,vmice)/100. ! to m/s ELSE ! By default, fallspeed velocity of ice crystals according to Heymsfield & Donner 1990 velo(i) = fallv_tun*cice_velo*((iwcg/1000.)**dice_velo) ENDIF ENDDO 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 :: 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 ! 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) .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. * 2./3. * tke(i) / 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 ) * 2./3. * tke(i) * 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) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! Calculate qsat following ECMWF method !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ USE yoethf_mod_h USE yomcst_mod_h IMPLICIT NONE include "FCTTRE.h" INTEGER, INTENT(IN) :: klon ! number of horizontal grid points REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature in K REAL, INTENT(IN), DIMENSION(klon) :: qtot ! total specific water in kg/kg REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa REAL, INTENT(IN) :: tref ! reference temperature in K LOGICAL, INTENT(IN) :: flagth ! flag for qsat calculation for thermals INTEGER, INTENT(IN) :: phase ! phase: 0=depend on temperature sign (temp>tref -> liquid, temp<tref, solid) ! 1=liquid ! 2=solid REAL, INTENT(OUT), DIMENSION(klon) :: qs ! saturation specific humidity [kg/kg] REAL, INTENT(OUT), DIMENSION(klon) :: dqs ! derivation of saturation specific humidity wrt T REAL delta, cor, cvm5 INTEGER i DO i=1,klon IF (phase .EQ. 1) THEN delta=0. ELSEIF (phase .EQ. 2) THEN delta=1. ELSE delta=MAX(0.,SIGN(1.,tref-temp(i))) ENDIF IF (flagth) THEN cvm5=R5LES*(1.-delta) + R5IES*delta ELSE cvm5 = R5LES*RLVTT*(1.-delta) + R5IES*RLSTT*delta cvm5 = cvm5 /RCPD/(1.0+RVTMP2*(qtot(i))) ENDIF qs(i)= R2ES*FOEEW(temp(i),delta)/pressure(i) qs(i)=MIN(0.5,qs(i)) cor=1./(1.-RETV*qs(i)) qs(i)=qs(i)*cor dqs(i)= FOEDE(temp(i),delta,cvm5,qs(i),cor) END DO END SUBROUTINE CALC_QSAT_ECMWF !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE CALC_GAMMASAT(klon,temp,qtot,pressure,gammasat,dgammasatdt) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! programme that calculates the gammasat parameter that determines the ! homogeneous condensation thresholds for cold (<0oC) clouds ! condensation at q>gammasat*qsat ! Etienne Vignon, March 2021 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ use lmdz_lscp_ini, only: iflag_gammasat, temp_nowater, RTT use lmdz_lscp_ini, only: a_homofreez, b_homofreez, delta_hetfreez IMPLICIT NONE INTEGER, INTENT(IN) :: klon ! number of horizontal grid points REAL, INTENT(IN), DIMENSION(klon) :: temp ! temperature in K REAL, INTENT(IN), DIMENSION(klon) :: qtot ! total specific water in kg/kg REAL, INTENT(IN), DIMENSION(klon) :: pressure ! pressure in Pa REAL, INTENT(OUT), DIMENSION(klon) :: gammasat ! coefficient to multiply qsat with to calculate saturation REAL, INTENT(OUT), DIMENSION(klon) :: dgammasatdt ! derivative of gammasat wrt temperature REAL, DIMENSION(klon) :: qsi,qsl,dqsl,dqsi REAL f_homofreez, fac INTEGER i CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,1,.false.,qsl,dqsl) CALL CALC_QSAT_ECMWF(klon,temp,qtot,pressure,RTT,2,.false.,qsi,dqsi) DO i = 1, klon IF ( temp(i) .GE. RTT ) THEN ! warm clouds: condensation at saturation wrt liquid gammasat(i) = 1. dgammasatdt(i) = 0. ELSE ! cold clouds: qsi > qsl ! homogeneous freezing of aerosols, according to ! Koop, 2000 and Ren and MacKenzie, 2005 (QJRMS) ! 'Cirrus regime' ! if f_homofreez > qsl / qsi, liquid nucleation ! if f_homofreez < qsl / qsi, homogeneous freezing of aerosols ! Note: f_homofreez = qsl / qsi for temp ~= -38degC f_homofreez = a_homofreez - temp(i) / b_homofreez IF ( iflag_gammasat .GE. 3 ) THEN ! condensation at homogeneous freezing threshold for temp < -38 degC ! condensation at liquid saturation for temp > -38 degC IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN gammasat(i) = f_homofreez dgammasatdt(i) = - 1. / b_homofreez ELSE gammasat(i) = qsl(i) / qsi(i) dgammasatdt(i) = ( dqsl(i) * qsi(i) - dqsi(i) * qsl(i) ) / qsi(i) / qsi(i) ENDIF ELSEIF ( iflag_gammasat .EQ. 2 ) THEN ! condensation at homogeneous freezing threshold for temp < -38 degC ! condensation at a threshold linearly decreasing between homogeneous ! freezing and ice saturation for -38 degC < temp < temp_nowater ! condensation at ice saturation for temp > temp_nowater ! If temp_nowater = 235.15 K, this is equivalent to iflag_gammasat = 1 IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN gammasat(i) = f_homofreez dgammasatdt(i) = - 1. / b_homofreez ELSEIF ( temp(i) .LE. temp_nowater ) THEN ! Here, we assume that f_homofreez = qsl / qsi for temp = -38 degC = 235.15 K dgammasatdt(i) = ( a_homofreez - 235.15 / b_homofreez - 1. ) & / ( 235.15 - temp_nowater ) gammasat(i) = dgammasatdt(i) * ( temp(i) - temp_nowater ) + 1. ELSE gammasat(i) = 1. dgammasatdt(i) = 0. ENDIF ELSEIF ( iflag_gammasat .EQ. 1 ) THEN ! condensation at homogeneous freezing threshold for temp < -38 degC ! condensation at ice saturation for temp > -38 degC IF ( f_homofreez .LE. qsl(i) / qsi(i) ) THEN gammasat(i) = f_homofreez dgammasatdt(i) = - 1. / b_homofreez ELSE gammasat(i) = 1. dgammasatdt(i) = 0. ENDIF ELSE ! condensation at ice saturation for temp < -38 degC ! condensation at ice saturation for temp > -38 degC gammasat(i) = 1. dgammasatdt(i) = 0. ENDIF ! Note that the delta_hetfreez parameter allows to linearly decrease the ! condensation threshold between the calculated threshold and the ice saturation ! for delta_hetfreez = 1, the threshold is the calculated condensation threshold ! for delta_hetfreez = 0, the threshold is the ice saturation gammasat(i) = ( 1. - delta_hetfreez ) + delta_hetfreez * gammasat(i) dgammasatdt(i) = delta_hetfreez * dgammasatdt(i) ENDIF END DO END SUBROUTINE CALC_GAMMASAT !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE DISTANCE_TO_CLOUD_TOP(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop) !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ USE lmdz_lscp_ini, ONLY : rd,rg,tresh_cl IMPLICIT NONE INTEGER, INTENT(IN) :: klon,klev !number of horizontal and vertical grid points INTEGER, INTENT(IN) :: k ! vertical index REAL, INTENT(IN), DIMENSION(klon,klev) :: temp ! temperature in K REAL, INTENT(IN), DIMENSION(klon,klev) :: pplay ! pressure middle layer in Pa REAL, INTENT(IN), DIMENSION(klon,klev+1) :: paprs ! pressure interfaces in Pa REAL, INTENT(IN), DIMENSION(klon,klev) :: rneb ! cloud fraction REAL, INTENT(OUT), DIMENSION(klon) :: distcltop1D ! distance from cloud top REAL, INTENT(OUT), DIMENSION(klon) :: temp_cltop ! temperature of cloud top REAL dzlay(klon,klev) REAL zlay(klon,klev) REAL dzinterf INTEGER i,k_top, kvert LOGICAL bool_cl DO i=1,klon ! Initialization height middle of first layer dzlay(i,1) = Rd * temp(i,1) / rg * log(paprs(i,1)/paprs(i,2)) zlay(i,1) = dzlay(i,1)/2 DO kvert=2,klev IF (kvert.EQ.klev) THEN dzlay(i,kvert) = 2*(rd * temp(i,kvert) / rg * log(paprs(i,kvert)/pplay(i,kvert))) ELSE dzlay(i,kvert) = rd * temp(i,kvert) / rg * log(paprs(i,kvert)/paprs(i,kvert+1)) ENDIF dzinterf = rd * temp(i,kvert) / rg * log(pplay(i,kvert-1)/pplay(i,kvert)) zlay(i,kvert) = zlay(i,kvert-1) + dzinterf ENDDO ENDDO DO i=1,klon k_top = k IF (rneb(i,k) .LE. tresh_cl) THEN bool_cl = .FALSE. ELSE bool_cl = .TRUE. ENDIF DO WHILE ((bool_cl) .AND. (k_top .LE. klev)) ! find cloud top IF (rneb(i,k_top) .GT. tresh_cl) THEN k_top = k_top + 1 ELSE bool_cl = .FALSE. k_top = k_top - 1 ENDIF ENDDO k_top=min(k_top,klev) !dist to top is dist between current layer and layer of cloud top (from middle to middle) + dist middle to !interf for layer of cloud top distcltop1D(i) = zlay(i,k_top) - zlay(i,k) + dzlay(i,k_top)/2 temp_cltop(i) = temp(i,k_top) ENDDO ! klon END SUBROUTINE DISTANCE_TO_CLOUD_TOP !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ FUNCTION GAMMAINC ( p, x ) !*****************************************************************************80 ! !! GAMMAINC computes the regularized lower incomplete Gamma Integral ! ! Modified: ! ! 20 January 2008 ! ! Author: ! ! Original FORTRAN77 version by B Shea. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! B Shea, ! Algorithm AS 239: ! Chi-squared and Incomplete Gamma Integral, ! Applied Statistics, ! Volume 37, Number 3, 1988, pages 466-473. ! ! Parameters: ! ! Input, real X, P, the parameters of the incomplete ! gamma ratio. 0 <= X, and 0 < P. ! ! Output, real GAMMAINC, the value of the incomplete ! Gamma integral. ! IMPLICIT NONE REAL A REAL AN REAL ARG REAL B REAL C REAL, PARAMETER :: ELIMIT = - 88.0E+00 REAL GAMMAINC REAL, PARAMETER :: OFLO = 1.0E+37 REAL P REAL, PARAMETER :: PLIMIT = 1000.0E+00 REAL PN1 REAL PN2 REAL PN3 REAL PN4 REAL PN5 REAL PN6 REAL RN REAL, PARAMETER :: TOL = 1.0E-14 REAL X REAL, PARAMETER :: XBIG = 1.0E+08 GAMMAINC = 0.0E+00 IF ( X == 0.0E+00 ) THEN GAMMAINC = 0.0E+00 RETURN END IF ! ! IF P IS LARGE, USE A NORMAL APPROXIMATION. ! IF ( PLIMIT < P ) THEN PN1 = 3.0E+00 * SQRT ( P ) * ( ( X / P )**( 1.0E+00 / 3.0E+00 ) & + 1.0E+00 / ( 9.0E+00 * P ) - 1.0E+00 ) GAMMAINC = 0.5E+00 * ( 1. + ERF ( PN1 ) ) RETURN END IF ! ! IF X IS LARGE SET GAMMAD = 1. ! IF ( XBIG < X ) THEN GAMMAINC = 1.0E+00 RETURN END IF ! ! USE PEARSON'S SERIES EXPANSION. ! (NOTE THAT P IS NOT LARGE ENOUGH TO FORCE OVERFLOW IN ALOGAM). ! IF ( X <= 1.0E+00 .OR. X < P ) THEN ARG = P * LOG ( X ) - X - LOG_GAMMA ( P + 1.0E+00 ) C = 1.0E+00 GAMMAINC = 1.0E+00 A = P DO A = A + 1.0E+00 C = C * X / A GAMMAINC = GAMMAINC + C IF ( C <= TOL ) THEN EXIT END IF END DO ARG = ARG + LOG ( GAMMAINC ) IF ( ELIMIT <= ARG ) THEN GAMMAINC = EXP ( ARG ) ELSE GAMMAINC = 0.0E+00 END IF ! ! USE A CONTINUED FRACTION EXPANSION. ! ELSE ARG = P * LOG ( X ) - X - LOG_GAMMA ( P ) A = 1.0E+00 - P B = A + X + 1.0E+00 C = 0.0E+00 PN1 = 1.0E+00 PN2 = X PN3 = X + 1.0E+00 PN4 = X * B GAMMAINC = PN3 / PN4 DO A = A + 1.0E+00 B = B + 2.0E+00 C = C + 1.0E+00 AN = A * C PN5 = B * PN3 - AN * PN1 PN6 = B * PN4 - AN * PN2 IF ( PN6 /= 0.0E+00 ) THEN RN = PN5 / PN6 IF ( ABS ( GAMMAINC - RN ) <= MIN ( TOL, TOL * RN ) ) THEN EXIT END IF GAMMAINC = RN END IF PN1 = PN3 PN2 = PN4 PN3 = PN5 PN4 = PN6 ! ! RE-SCALE TERMS IN CONTINUED FRACTION IF TERMS ARE LARGE. ! IF ( OFLO <= ABS ( PN5 ) ) THEN PN1 = PN1 / OFLO PN2 = PN2 / OFLO PN3 = PN3 / OFLO PN4 = PN4 / OFLO END IF END DO ARG = ARG + LOG ( GAMMAINC ) IF ( ELIMIT <= ARG ) THEN GAMMAINC = 1.0E+00 - EXP ( ARG ) ELSE GAMMAINC = 1.0E+00 END IF END IF RETURN END FUNCTION GAMMAINC !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ END MODULE lmdz_lscp_tools