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 == 1) THEN ! so-called 'empirical parameterization' in Stubenrauch et al. 2019 IF (tempc >= -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 == 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 lmdz_print_control, 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 USE lmdz_abort_physic, ONLY: abort_physic 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<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 == 1).OR.(iflag_ice_thermo >= 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 == 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) >0.) THEN dicefracdT(i) = exposant_glace * (icefrac(i)**(exposant_glace - 1.)) & / (t_glace_min - t_glace_max) ENDIF IF ((icefrac(i)==0).OR.(icefrac(i)==1)) THEN dicefracdT(i) = 0. ENDIF ENDIF ! function of temperature used in CMIP6 physics IF (iflag_t_glace == 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) >0.) .AND. (liqfrac_tmp > 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 == 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 <=0) .OR. (liqfrac_tmp >= 1)) THEN dicefracdT(i) = 0. ENDIF ENDIF ! with CMIP6 function of temperature at cloud top IF ((iflag_t_glace == 5) .OR. (iflag_t_glace == 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 <=0) .OR. (liqfrac_tmp >= 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 == 6) THEN IF (temp(i) > 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 <=0) .OR. (liqfrac_tmp >= 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 of cloud top <-40°C, IF (iflag_t_glace >= 4) THEN IF ((temp_cltop(i) <= temp_nowater) .AND. (temp(i) <= t_glace_max)) THEN icefrac(i) = 1. dicefracdT(i) = 0. ENDIF ENDIF ENDDO ! klon END SUBROUTINE ICEFRAC_LSCP !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE ICEFRAC_LSCP_TURB(klon, dtime, temp, pplay, paprsdn, paprsup, qice_ini, snowcld, qtot_incl, cldfra, tke, tke_dissip, 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 (30/08/24) !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 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: tau_mixenv, lmix_mpc, naero5, gamma_snwretro, gamma_taud, capa_crystal USE lmdz_lscp_ini, ONLY: eps IMPLICIT NONE INTEGER, INTENT(IN) :: klon !--number of horizontal grid points REAL, INTENT(IN) :: dtime !--time step [s] 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) :: qtot_incl !--specific total cloud water content, 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 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(OUT), DIMENSION(klon) :: icefrac !--fraction of ice in condensed water [-] REAL, INTENT(OUT), DIMENSION(klon) :: dicefracdT REAL, INTENT(OUT), DIMENSION(klon) :: cldfraliq !--fraction of cldfra which is liquid only REAL, INTENT(OUT), DIMENSION(klon) :: sigma2_icefracturb !--Temporary REAL, INTENT(OUT), DIMENSION(klon) :: mean_icefracturb !--Temporary 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 :: qsnowcld_incl !REAL :: capa_crystal !--Capacitance of ice crystals [-] REAL :: water_vapor_diff !--Water-vapour diffusion coefficient in air [m2/s] (function of T&P) REAL :: air_thermal_conduct !--Thermal conductivity of air [J/m/K/s] (function of T) REAL :: C0 !--Lagrangian structure function [-] REAL :: tau_mixingenv REAL :: tau_dissipturb REAL :: invtau_phaserelax REAL :: sigma2_pdf, mean_pdf REAL :: ai, bi, B0 REAL :: sursat_iceliq REAL :: sursat_env 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 :: rho_ice !--ice density [kg/m3] REAL :: cldfra1D REAL :: deltaz, rho_air REAL :: psati !--saturation vapor pressure wrt i [Pa] C0 = 10. !--value assumed in Field2014 rho_ice = 950. sursat_iceext = -0.1 !capa_crystal = 1. !r_ice qzero(:) = 0. cldfraliq(:) = 0. icefrac(:) = 0. dicefracdT(:) = 0. sigma2_icefracturb(:) = 0. mean_icefracturb(:) = 0. !--wrt liquid water 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 !deltaz = ( paprsdn(i) - paprsup(i) ) / RG / rho_air(i) ! because cldfra is intent in, but can be locally modified due to test cldfra1D = cldfra(i) IF (cldfra(i) <= 0.) THEN qvap_cld(i) = 0. qliq(i) = 0. qice(i) = 0. cldfraliq(i) = 0. icefrac(i) = 0. dicefracdT(i) = 0. ! If there is a cloud ELSE IF (cldfra(i) >= 1.0) THEN cldfra1D = 1.0 END IF ! T>0°C, no ice allowed IF (temp(i) >= 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) <= 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. ! MPC temperature ELSE ! Not enough TKE IF (tke_dissip(i) <= eps) THEN 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. ! Enough TKE ELSE print*,"MOUCHOIRACTIVE" !--------------------------------------------------------- !-- ICE SUPERSATURATION PDF !--------------------------------------------------------- !--If -38°C< T <0°C and there is enough turbulence, !--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. sursat_iceliq = qsatl(i) / qsati(i) - 1. psati = qsati(i) * pplay(i) / (RD / RV) !-------------- MICROPHYSICAL TERMS -------------- !--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 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 qiceini_incl = qice_ini(i) / cldfra1D qsnowcld_incl = snowcld(i) * RG * dtime / (paprsdn(i) - paprsup(i)) / cldfra1D sursat_env = max(0., (qtot_incl(i) - qiceini_incl) / qsati(i) - 1.) IF (qiceini_incl > eps) THEN nb_crystals = 1.e3 * 5.94e-5 * (RTT - temp(i))**3.33 * naero5**(0.0264 * (RTT - temp(i)) + 0.0033) lambda_PSD = ((RPI * rho_ice * nb_crystals ) / (rho_air * (qiceini_incl + gamma_snwretro * qsnowcld_incl))) ** (1. / 3.) N0_PSD = nb_crystals * lambda_PSD moment1_PSD = N0_PSD / lambda_PSD**2 ELSE moment1_PSD = 0. ENDIF !--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) ! Old way of estimating moment1 : spherical crystals + monodisperse PSD ! nb_crystals = rho_air * qiceini_incl / ( 4. / 3. * RPI * r_ice**3. * rho_ice ) ! moment1_PSD = nb_crystals * r_ice !----------------- TURBULENT SOURCE/SINK TERMS ----------------- !--Tau_mixingenv is the time needed to homogeneize the parcel !--with its environment by turbulent diffusion over the parcel !--length scale !--if lmix_mpc <0, tau_mixigenv value is prescribed !--else tau_mixigenv value is derived from tke_dissip and lmix_mpc !--Tau_dissipturb is the time needed turbulence to decay due to !--viscosity ai = RG / RD / temp(i) * (RD * RLSTT / RCPD / RV / temp(i) - 1.) IF (lmix_mpc > 0) THEN tau_mixingenv = (lmix_mpc**2. / tke_dissip(i))**(1. / 3.) ELSE tau_mixingenv = tau_mixenv ENDIF tau_dissipturb = gamma_taud * 2. * 2. / 3. * tke(i) / tke_dissip(i) / C0 !--------------------- PDF COMPUTATIONS --------------------- !--Formulae for sigma2_pdf (variance), mean of PDF in Furtado2016 !--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 liqfra_max = MAX(0., (MIN (1., (qtot_incl(i) - qiceini_incl - qsati(i) * (1 + sursat_iceext)) / (qsatl(i) - qsati(i))))) sigma2_pdf = 1. / 2. * (ai**2) * 2. / 3. * tke(i) * tau_dissipturb / (invtau_phaserelax + 1. / tau_mixingenv) mean_pdf = sursat_env * 1. / tau_mixingenv / (invtau_phaserelax + 1. / tau_mixingenv) cldfraliq(i) = 0.5 * (1. - erf((sursat_iceliq - mean_pdf) / (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 - mean_pdf)**2. / (2. * sigma2_pdf)) & - qsati(i) * cldfraliq(i) * (sursat_iceliq - mean_pdf) sigma2_icefracturb(i) = sigma2_pdf mean_icefracturb(i) = mean_pdf !------------ SPECIFIC VAPOR CONTENT AND WATER CONSERVATION ------------ IF ((qliq_incl <= eps) .OR. (cldfraliq(i) <= 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 >= qtot_incl(i)) THEN qvap_incl = qsati(i) qliq_incl = qtot_incl(i) - qvap_incl qice_incl = 0. ELSEIF ((qvap_incl + qliq_incl) >= 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 icefrac(i) = qice(i) / (qice(i) + qliq(i)) dicefracdT(i) = 0. !PRINT*,'MPC turb' END IF ! Enough TKE END IF ! MPC temperature END IF ! 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 lmdz_yoethf USE lmdz_yomcst 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, tempgammasat*qsat ! Etienne Vignon, March 2021 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ USE lmdz_lscp_ini, ONLY: iflag_gammasat, t_glace_min, RTT 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 fcirrus, fac REAL, PARAMETER :: acirrus = 2.349 REAL, PARAMETER :: bcirrus = 259.0 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) >= RTT) THEN ! warm clouds: condensation at saturation wrt liquid gammasat(i) = 1. dgammasatdt(i) = 0. ELSEIF ((temp(i) < RTT) .AND. (temp(i) > t_glace_min)) THEN IF (iflag_gammasat >= 2) THEN gammasat(i) = qsl(i) / qsi(i) dgammasatdt(i) = (dqsl(i) * qsi(i) - dqsi(i) * qsl(i)) / qsi(i) / qsi(i) ELSE gammasat(i) = 1. dgammasatdt(i) = 0. ENDIF ELSE IF (iflag_gammasat >=1) THEN ! homogeneous freezing of aerosols, according to ! Koop, 2000 and Karcher 2008, QJRMS ! 'Cirrus regime' fcirrus = acirrus - temp(i) / bcirrus IF (fcirrus > qsl(i) / qsi(i)) THEN gammasat(i) = qsl(i) / qsi(i) dgammasatdt(i) = (dqsl(i) * qsi(i) - dqsi(i) * qsl(i)) / qsi(i) / qsi(i) ELSE gammasat(i) = fcirrus dgammasatdt(i) = -1.0 / bcirrus ENDIF ELSE gammasat(i) = 1. dgammasatdt(i) = 0. ENDIF 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==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) <= tresh_cl) THEN bool_cl = .FALSE. ELSE bool_cl = .TRUE. ENDIF DO WHILE ((bool_cl) .AND. (k_top <= klev)) ! find cloud top IF (rneb(i, k_top) > 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 !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ END MODULE lmdz_lscp_tools