MODULE lmdz_lscp_precip !---------------------------------------------------------------- ! Module for the process-oriented treament of precipitation ! that are called in LSCP ! Authors: Atelier Nuage (G. Riviere, L. Raillard, M. Wimmer, ! N. Dutrievoz, E. Vignon, A. Borella, et al.) ! Jan. 2024 IMPLICIT NONE CONTAINS !---------------------------------------------------------------- ! historical (till CMIP6) version of the pre-cloud formation ! precipitation scheme containing precip evaporation and melting SUBROUTINE histprecip_precld( & klon, dtime, iftop, paprsdn, paprsup, pplay, zt, ztupnew, zq, & zmqc, zneb, znebprecip, znebprecipclr, & zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, dqreva, dqssub & ) USE lmdz_lscp_ini, ONLY : iflag_evap_prec USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub, ztfondue USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG 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) :: iftop !--if top of the column 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) :: pplay !--pressure in the middle of the layer [Pa] REAL, INTENT(INOUT), DIMENSION(klon) :: zt !--current temperature [K] REAL, INTENT(IN), DIMENSION(klon) :: ztupnew !--updated temperature of the overlying layer [K] REAL, INTENT(INOUT), DIMENSION(klon) :: zq !--current water vapor specific humidity (includes evaporated qi and ql) [kg/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: zmqc !--specific humidity in the precipitation falling from the upper layer [kg/kg] REAL, INTENT(IN), DIMENSION(klon) :: zneb !--cloud fraction IN THE LAYER ABOVE [-] REAL, INTENT(INOUT), DIMENSION(klon) :: znebprecip !--fraction of precipitation IN THE LAYER ABOVE [-] REAL, INTENT(INOUT), DIMENSION(klon) :: znebprecipclr !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-] REAL, INTENT(INOUT), DIMENSION(klon) :: zrfl !--flux of rain gridbox-mean coming from the layer above [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: zrflclr !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: zrflcld !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: zifl !--flux of snow gridbox-mean coming from the layer above [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: ziflclr !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: ziflcld !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2] REAL, INTENT(OUT), DIMENSION(klon) :: dqreva !--rain tendency due to evaporation [kg/kg/s] REAL, INTENT(OUT), DIMENSION(klon) :: dqssub !--snow tendency due to sublimation [kg/kg/s] REAL :: zmair REAL :: zcpair, zcpeau REAL, DIMENSION(klon) :: qzero REAL, DIMENSION(klon) :: zqs, zdqs REAL, DIMENSION(klon) :: qsl, qsi REAL, DIMENSION(klon) :: dqsl, dqsi REAL :: zqev0 REAL :: zqev, zqevt REAL :: zqevi, zqevti REAL, DIMENSION(klon) :: zrfln, zifln REAL :: zmelt INTEGER :: i qzero(:) = 0. ! -------------------------------------------------------------------- ! P1> Thermalization of precipitation falling from the overlying layer ! -------------------------------------------------------------------- ! Computes air temperature variation due to enthalpy transported by ! precipitation. Precipitation is then thermalized with the air in the ! layer. ! The precipitation should remain thermalized throughout the different ! thermodynamical transformations. ! The corresponding water mass should ! be added when calculating the layer's enthalpy change with ! temperature ! --------------------------------------------------------------------- IF (iftop) THEN DO i = 1, klon zmqc(i) = 0. ENDDO ELSE DO i = 1, klon zmair=(paprsdn(i)-paprsup(i))/RG ! no condensed water so cp=cp(vapor+dry air) ! RVTMP2=rcpv/rcpd-1 zcpair=RCPD*(1.0+RVTMP2*zq(i)) zcpeau=RCPD*RVTMP2 ! zmqc: precipitation mass that has to be thermalized with ! layer's air so that precipitation at the ground has the ! same temperature as the lowermost layer zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer zt(i) = ( ztupnew(i)*zmqc(i)*zcpeau + zcpair*zt(i) ) & / (zcpair + zmqc(i)*zcpeau) ENDDO ENDIF ! -------------------------------------------------------------------- ! P2> Precipitation evaporation/sublimation/melting ! -------------------------------------------------------------------- ! A part of the precipitation coming from above is evaporated/sublimated/melted. ! -------------------------------------------------------------------- IF (iflag_evap_prec.GE.1) THEN ! Calculation of saturation specific humidity ! depending on temperature: CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:),RTT,0,.false.,zqs,zdqs) ! wrt liquid water CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:),RTT,1,.false.,qsl,dqsl) ! wrt ice CALL calc_qsat_ecmwf(klon,zt,qzero,pplay(:),RTT,2,.false.,qsi,dqsi) DO i = 1, klon ! if precipitation IF (zrfl(i)+zifl(i).GT.0.) THEN ! LudoTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec>=4). ! c_iso: likely important to distinguish cs from neb isotope precipitation IF (iflag_evap_prec.GE.4) THEN zrfl(i) = zrflclr(i) zifl(i) = ziflclr(i) ENDIF IF (iflag_evap_prec.EQ.1) THEN znebprecip(i)=zneb(i) ELSE znebprecip(i)=MAX(zneb(i),znebprecip(i)) ENDIF IF (iflag_evap_prec.GT.4) THEN ! Max evaporation not to saturate the clear sky precip fraction ! i.e. the fraction where evaporation occurs zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecipclr(i)) ELSEIF (iflag_evap_prec .EQ. 4) THEN ! Max evaporation not to saturate the whole mesh ! Pay attention -> lead to unrealistic and excessive evaporation zqev0 = MAX(0.0, zqs(i)-zq(i)) ELSE ! Max evap not to saturate the fraction below the cloud zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecip(i)) ENDIF ! Evaporation of liquid precipitation coming from above ! dP/dz=beta*(1-q/qsat)*sqrt(P) ! formula from Sundquist 1988, Klemp & Wilhemson 1978 ! LTP: evaporation only in the clear sky part (iflag_evap_prec>=4) IF (iflag_evap_prec.EQ.3) THEN zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl(i)) & *SQRT(zrfl(i)/max(1.e-4,znebprecip(i))) & *(paprsdn(i)-paprsup(i))/pplay(i)*zt(i)*RD/RG ELSE IF (iflag_evap_prec.GE.4) THEN zqevt = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsl(i)) & *SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) & *(paprsdn(i)-paprsup(i))/pplay(i)*zt(i)*RD/RG ELSE zqevt = 1.*coef_eva*(1.0-zq(i)/qsl(i))*SQRT(zrfl(i)) & *(paprsdn(i)-paprsup(i))/pplay(i)*zt(i)*RD/RG ENDIF zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) * RG*dtime/(paprsdn(i)-paprsup(i)) ! sublimation of the solid precipitation coming from above IF (iflag_evap_prec.EQ.3) THEN zqevti = znebprecip(i)*coef_sub*(1.0-zq(i)/qsi(i)) & *SQRT(zifl(i)/max(1.e-4,znebprecip(i))) & *(paprsdn(i)-paprsup(i))/pplay(i)*zt(i)*RD/RG ELSE IF (iflag_evap_prec.GE.4) THEN zqevti = znebprecipclr(i)*coef_sub*(1.0-zq(i)/qsi(i)) & *SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) & *(paprsdn(i)-paprsup(i))/pplay(i)*zt(i)*RD/RG ELSE zqevti = 1.*coef_sub*(1.0-zq(i)/qsi(i))*SQRT(zifl(i)) & *(paprsdn(i)-paprsup(i))/pplay(i)*zt(i)*RD/RG ENDIF zqevti = MAX(0.0,MIN(zqevti,zifl(i))) * RG*dtime/(paprsdn(i)-paprsup(i)) ! A. JAM ! Evaporation limit: we ensure that the layer's fraction below ! the cloud or the whole mesh (depending on iflag_evap_prec) ! does not reach saturation. In this case, we ! redistribute zqev0 conserving the ratio liquid/ice IF (zqevt+zqevti.GT.zqev0) THEN zqev=zqev0*zqevt/(zqevt+zqevti) zqevi=zqev0*zqevti/(zqevt+zqevti) ELSE zqev=zqevt zqevi=zqevti ENDIF !--Diagnostics dqreva(i) = - zqev / dtime dqssub(i) = - zqevti / dtime ! New solid and liquid precipitation fluxes after evap and sublimation zrfln(i) = Max(0.,zrfl(i) - zqev*(paprsdn(i)-paprsup(i))/RG/dtime) zifln(i) = Max(0.,zifl(i) - zqevi*(paprsdn(i)-paprsup(i))/RG/dtime) ! vapor, temperature, precip fluxes update ! vapor is updated after evaporation/sublimation (it is increased) zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) & * (RG/(paprsdn(i)-paprsup(i)))*dtime ! zmqc is the total condensed water in the precip flux (it is decreased) zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) & * (RG/(paprsdn(i)-paprsup(i)))*dtime ! air and precip temperature (i.e., gridbox temperature) ! is updated due to latent heat cooling zt(i) = zt(i) + (zrfln(i)-zrfl(i)) & * (RG/(paprsdn(i)-paprsup(i)))*dtime & * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) & + (zifln(i)-zifl(i)) & * (RG/(paprsdn(i)-paprsup(i)))*dtime & * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) ! New values of liquid and solid precipitation zrfl(i) = zrfln(i) zifl(i) = zifln(i) ! c_iso here call_reevap that updates isotopic zrfl, zifl (in inout) ! due to evap + sublim IF (iflag_evap_prec.GE.4) THEN zrflclr(i) = zrfl(i) ziflclr(i) = zifl(i) IF(zrflclr(i) + ziflclr(i).LE.0) THEN znebprecipclr(i) = 0.0 ENDIF zrfl(i) = zrflclr(i) + zrflcld(i) zifl(i) = ziflclr(i) + ziflcld(i) ENDIF ! c_iso duplicate for isotopes or loop on isotopes ! Melting: zmelt = ((zt(i)-RTT)/(ztfondue-RTT)) ! JYG ! precip fraction that is melted zmelt = MIN(MAX(zmelt,0.),1.) ! update of rainfall and snowfall due to melting IF (iflag_evap_prec.GE.4) THEN zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i) zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i) zrfl(i)=zrflclr(i)+zrflcld(i) ELSE zrfl(i)=zrfl(i)+zmelt*zifl(i) ENDIF ! c_iso: melting of isotopic precipi with zmelt (no fractionation) ! Latent heat of melting because of precipitation melting ! NB: the air + precip temperature is simultaneously updated zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprsdn(i)-paprsup(i)) & *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) IF (iflag_evap_prec.GE.4) THEN ziflclr(i)=ziflclr(i)*(1.-zmelt) ziflcld(i)=ziflcld(i)*(1.-zmelt) zifl(i)=ziflclr(i)+ziflcld(i) ELSE zifl(i)=zifl(i)*(1.-zmelt) ENDIF ELSE ! if no precip, we reinitialize the cloud fraction used for the precip to 0 znebprecip(i)=0. ENDIF ! (zrfl(i)+zifl(i).GT.0.) ENDDO ! loop on klon ENDIF ! (iflag_evap_prec>=1) END SUBROUTINE histprecip_precld SUBROUTINE histprecip_postcld( & klon, dtime, iftop, paprsdn, paprsup, pplay, ctot_vol, ptconv, zdqsdT_raw, & zt, zq, zoliq, zoliql, zoliqi, zcond, zfice, zmqc, & rneb, znebprecipclr, znebprecipcld, zneb, tot_zneb, zrho_up, zvelo_up, & zrfl, zrflclr, zrflcld, zifl, ziflclr, ziflcld, & zradocond, zradoice, dqrauto, dqsauto & ) USE lmdz_lscp_ini, ONLY : RD, RG, RCPD, RVTMP2, RLSTT, RLMLT, RTT USE lmdz_lscp_ini, ONLY : cld_lc_con, cld_tau_con, cld_expo_con, ffallv_con, & cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, ffallv_lsc, & seuil_neb, rain_int_min, cice_velo, dice_velo USE lmdz_lscp_ini, ONLY : iflag_evap_prec, iflag_cloudth_vert, iflag_rain_incloud_vol, & iflag_autoconversion, ok_radocond_snow, ok_bug_phase_lscp, & niter_lscp USE lmdz_lscp_tools, ONLY : fallice_velocity IMPLICIT NONE INTEGER, INTENT(IN) :: klon !--number of horizontal grid points [-] REAL, INTENT(IN) :: dtime !--time step [s] LOGICAL, INTENT(IN) :: iftop !--if top of the column 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) :: pplay !--pressure in the middle of the layer [Pa] REAL, INTENT(IN), DIMENSION(klon) :: ctot_vol !--volumic cloud fraction [-] LOGICAL, INTENT(IN), DIMENSION(klon) :: ptconv !--true if we are in a convective point REAL, INTENT(IN), DIMENSION(klon) :: zdqsdT_raw !--derivative of qsat w.r.t. temperature times L/Cp [SI] REAL, INTENT(INOUT), DIMENSION(klon) :: zt !--current temperature [K] REAL, INTENT(INOUT), DIMENSION(klon) :: zq !--current water vapor specific humidity [kg/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: zoliq !--current liquid+ice water specific humidity [kg/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: zoliql !--current liquid water specific humidity [kg/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: zoliqi !--current ice water specific humidity [kg/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: zcond !--liquid+ice water specific humidity AFTER CONDENSATION (no precip process) [kg/kg] REAL, INTENT(IN), DIMENSION(klon) :: zfice !--ice fraction AFTER CONDENSATION [-] REAL, INTENT(IN), DIMENSION(klon) :: zmqc !--specific humidity in the precipitation falling from the upper layer [kg/kg] REAL, INTENT(IN), DIMENSION(klon) :: rneb !--cloud fraction [-] REAL, INTENT(INOUT), DIMENSION(klon) :: znebprecipclr !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-] REAL, INTENT(INOUT), DIMENSION(klon) :: znebprecipcld !--fraction of precipitation in the cloud IN THE LAYER ABOVE [-] REAL, INTENT(INOUT), DIMENSION(klon) :: zneb !--cloud fraction IN THE LAYER ABOVE [-] REAL, INTENT(INOUT), DIMENSION(klon) :: tot_zneb !--total cloud cover above the current mesh [-] REAL, INTENT(INOUT), DIMENSION(klon) :: zrho_up !--air density IN THE LAYER ABOVE [kg/m3] REAL, INTENT(INOUT), DIMENSION(klon) :: zvelo_up !--ice fallspeed velocity IN THE LAYER ABOVE [m/s] REAL, INTENT(INOUT), DIMENSION(klon) :: zrfl !--flux of rain gridbox-mean [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: zrflclr !--flux of rain gridbox-mean in clear sky [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: zrflcld !--flux of rain gridbox-mean in cloudy air [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: zifl !--flux of snow gridbox-mean [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: ziflclr !--flux of snow gridbox-mean in clear sky [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: ziflcld !--flux of snow gridbox-mean in cloudy air [kg/s/m2] REAL, INTENT(OUT), DIMENSION(klon) :: zradocond !--condensed water used in the radiation scheme [kg/kg] REAL, INTENT(OUT), DIMENSION(klon) :: zradoice !--condensed ice water used in the radiation scheme [kg/kg] REAL, INTENT(OUT), DIMENSION(klon) :: dqrauto !--rain tendency due to autoconversion of cloud liquid [kg/kg/s] REAL, INTENT(OUT), DIMENSION(klon) :: dqsauto !--snow tendency due to autoconversion of cloud ice [kg/kg/s] ! Local variables for precip fraction update REAL :: smallestreal REAL, DIMENSION(klon) :: tot_znebn, d_tot_zneb REAL, DIMENSION(klon) :: d_znebprecip_cld_clr, d_znebprecip_clr_cld REAL, DIMENSION(klon) :: d_zrfl_cld_clr, d_zifl_cld_clr REAL, DIMENSION(klon) :: d_zrfl_clr_cld, d_zifl_clr_cld ! Local variables for autoconversion REAL :: zct, zcl, zexpo, ffallv REAL :: zchau, zfroi REAL :: zrain, zsnow, zprecip REAL :: effective_zneb REAL, DIMENSION(klon) :: zrho, zvelo REAL, DIMENSION(klon) :: zdz, iwc ! Local variables for Bergeron process REAL :: zcp, coef1, DeltaT, Deltaq, Deltaqprecl REAL, DIMENSION(klon) :: zqpreci, zqprecl ! Local variables for calculation of quantity of condensates seen by the radiative scheme REAL, DIMENSION(klon) :: zradoliq REAL, DIMENSION(klon) :: ziflprev ! Misc INTEGER :: i, n ! Initialisation smallestreal=1.e-9 zqprecl(:) = 0. zqpreci(:) = 0. ziflprev(:) = 0. IF (iflag_evap_prec .GE. 4) THEN !Partitionning between precipitation coming from clouds and that coming from CS !0) Calculate tot_zneb, total cloud fraction above the cloud !assuming a maximum-random overlap (voir Jakob and Klein, 2000) DO i=1, klon tot_znebn(i) = 1. - (1.-tot_zneb(i))*(1 - max(rneb(i),zneb(i))) & /(1.-min(zneb(i),1.-smallestreal)) d_tot_zneb(i) = tot_znebn(i) - tot_zneb(i) tot_zneb(i) = tot_znebn(i) !1) Cloudy to clear air d_znebprecip_cld_clr(i) = znebprecipcld(i) - min(rneb(i),znebprecipcld(i)) IF (znebprecipcld(i) .GT. 0.) THEN d_zrfl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*zrflcld(i) d_zifl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*ziflcld(i) ELSE d_zrfl_cld_clr(i) = 0. d_zifl_cld_clr(i) = 0. ENDIF !2) Clear to cloudy air d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i) - d_tot_zneb(i) - zneb(i))) IF (znebprecipclr(i) .GT. 0) THEN d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i) d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*ziflclr(i) ELSE d_zrfl_clr_cld(i) = 0. d_zifl_clr_cld(i) = 0. ENDIF !Update variables znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i) znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i) zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i) ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i) zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i) ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i) ! c_iso do the same thing for isotopes precip ENDDO ENDIF ! Autoconversion ! ------------------------------------------------------------------------------- ! Initialisation DO i = 1, klon zrho(i) = pplay(i) / zt(i) / RD zdz(i) = (paprsdn(i)-paprsup(i)) / (zrho(i)*RG) iwc(i) = 0. zneb(i) = MAX(rneb(i),seuil_neb) ! variables for calculation of quantity of condensates seen by the radiative scheme ! NB. only zradocond and zradoice are outputed, but zradoliq is used if ok_radocond_snow ! is TRUE zradocond(i) = zoliq(i)/REAL(niter_lscp+1) zradoliq(i) = zoliq(i)*(1.-zfice(i))/REAL(niter_lscp+1) zradoice(i) = zoliq(i)*zfice(i)/REAL(niter_lscp+1) ENDDO DO n = 1, niter_lscp DO i=1,klon IF (rneb(i).GT.0.0) THEN iwc(i) = zrho(i) * zoliqi(i) / zneb(i) ! in-cloud ice condensate content ENDIF ENDDO CALL fallice_velocity(klon, iwc, zt, zrho, pplay, ptconv, zvelo) DO i = 1, klon IF (rneb(i).GT.0.0) THEN ! Initialization of zrain, zsnow and zprecip: zrain=0. zsnow=0. zprecip=0. ! c_iso same init for isotopes. Externalisation? IF (zneb(i).EQ.seuil_neb) THEN zprecip = 0.0 zsnow = 0.0 zrain= 0.0 ELSE IF (ptconv(i)) THEN ! if convective point zcl=cld_lc_con zct=1./cld_tau_con zexpo=cld_expo_con ffallv=ffallv_con ELSE zcl=cld_lc_lsc zct=1./cld_tau_lsc zexpo=cld_expo_lsc ffallv=ffallv_lsc ENDIF ! if vertical heterogeneity is taken into account, we use ! the "true" volume fraction instead of a modified ! surface fraction (which is larger and artificially ! reduces the in-cloud water). IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN effective_zneb=ctot_vol(i) ELSE effective_zneb=zneb(i) ENDIF ! Liquid water quantity to remove: zchau (Sundqvist, 1978) ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2) !......................................................... IF (iflag_autoconversion .EQ. 2) THEN ! two-steps resolution with niter_lscp=1 sufficient ! we first treat the second term (with exponential) in an explicit way ! and then treat the first term (-q/tau) in an exact way zchau=zoliql(i)*(1.-exp(-dtime/REAL(niter_lscp)*zct & *(1.-exp(-(zoliql(i)/effective_zneb/zcl)**zexpo)))) ELSE ! old explicit resolution with subtimesteps zchau = zct*dtime/REAL(niter_lscp)*zoliql(i) & *(1.0-EXP(-(zoliql(i)/effective_zneb/zcl)**zexpo)) ENDIF ! Ice water quantity to remove (Zender & Kiehl, 1997) ! dqice/dt=1/rho*d(rho*wice*qice)/dz !.................................... IF (iflag_autoconversion .EQ. 2) THEN ! exact resolution, niter_lscp=1 is sufficient but works only ! with iflag_vice=0 IF (zoliqi(i) .GT. 0.) THEN zfroi=(zoliqi(i)-((zoliqi(i)**(-dice_velo)) & +dice_velo*dtime/REAL(niter_lscp)*cice_velo/zdz(i)/zneb(i)*zrho(i)*ffallv)**(-1./dice_velo)) ELSE zfroi=0. ENDIF ELSE ! old explicit resolution with subtimesteps zfroi = dtime/REAL(niter_lscp)/zdz(i)*zoliqi(i)*zvelo(i) ENDIF zrain = MIN(MAX(zchau,0.0),zoliql(i)) zsnow = MIN(MAX(zfroi,0.0),zoliqi(i)) zprecip = MAX(zrain + zsnow,0.0) ENDIF IF (iflag_autoconversion .GE. 1) THEN ! debugged version with phase conservation through the autoconversion process zoliql(i) = MAX(zoliql(i)-1.*zrain , 0.0) zoliqi(i) = MAX(zoliqi(i)-1.*zsnow , 0.0) zoliq(i) = MAX(zoliq(i)-zprecip , 0.0) ELSE ! bugged version with phase resetting zoliql(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zrain , 0.0) zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zsnow , 0.0) zoliq(i) = MAX(zoliq(i)-zprecip , 0.0) ENDIF ! c_iso: call isotope_conversion (for readibility) or duplicate ! variables for calculation of quantity of condensates seen by the radiative scheme zradocond(i) = zradocond(i) + zoliq(i)/REAL(niter_lscp+1) zradoliq(i) = zradoliq(i) + zoliql(i)/REAL(niter_lscp+1) zradoice(i) = zradoice(i) + zoliqi(i)/REAL(niter_lscp+1) !--Diagnostics dqrauto(i) = dqrauto(i) + zrain / dtime dqsauto(i) = dqsauto(i) + zsnow / dtime ENDIF ! rneb >0 ENDDO ! i = 1,klon ENDDO ! n = 1,niter ! Precipitation flux calculation DO i = 1, klon IF (iflag_evap_prec.GE.4) THEN ziflprev(i)=ziflcld(i) ELSE ziflprev(i)=zifl(i)*zneb(i) ENDIF IF (rneb(i) .GT. 0.0) THEN ! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux: ! If T<0C, liquid precip are converted into ice, which leads to an increase in ! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly ! taken into account through a linearization of the equations and by approximating ! the liquid precipitation process with a "threshold" process. We assume that ! condensates are not modified during this operation. Liquid precipitation is ! removed (in the limit DeltaT<273.15-T). Solid precipitation increases. ! Water vapor increases as well ! Note that compared to fisrtilp, we always assume iflag_bergeron=2 IF (ok_bug_phase_lscp) THEN zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i) zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i)) ELSE zqpreci(i)=zcond(i)*zfice(i)-zoliqi(i) zqprecl(i)=zcond(i)*(1.-zfice(i))-zoliql(i) ENDIF zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) coef1 = rneb(i)*RLSTT/zcp*zdqsdT_raw(i) ! Computation of DT if all the liquid precip freezes DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1)) ! T should not exceed the freezing point ! that is Delta > RTT-zt(i) DeltaT = max( min( RTT-zt(i), DeltaT) , 0. ) zt(i) = zt(i) + DeltaT ! water vaporization due to temp. increase Deltaq = rneb(i)*zdqsdT_raw(i)*DeltaT ! we add this vaporized water to the total vapor and we remove it from the precip zq(i) = zq(i) + Deltaq ! The three "max" lines herebelow protect from rounding errors zcond(i) = max( zcond(i)- Deltaq, 0. ) ! liquid precipitation converted to ice precip Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. ) ! iced water budget zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.) ! c_iso : mv here condensation of isotopes + redispatchage en precip IF (iflag_evap_prec.GE.4) THEN zrflcld(i) = zrflcld(i)+zqprecl(i) & *(paprsdn(i)-paprsup(i))/(RG*dtime) ziflcld(i) = ziflcld(i)+ zqpreci(i) & *(paprsdn(i)-paprsup(i))/(RG*dtime) znebprecipcld(i) = rneb(i) zrfl(i) = zrflcld(i) + zrflclr(i) zifl(i) = ziflcld(i) + ziflclr(i) ELSE zrfl(i) = zrfl(i)+ zqprecl(i) & *(paprsdn(i)-paprsup(i))/(RG*dtime) zifl(i) = zifl(i)+ zqpreci(i) & *(paprsdn(i)-paprsup(i))/(RG*dtime) ENDIF ! c_iso : same for isotopes ENDIF ! rneb>0 ENDDO ! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min ! if iflag_evap_prec>=4 IF (iflag_evap_prec.GE.4) THEN DO i=1,klon IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN znebprecipclr(i) = min(znebprecipclr(i),max( & zrflclr(i)/ (MAX(znebprecipclr(i),seuil_neb)*rain_int_min), & ziflclr(i)/ (MAX(znebprecipclr(i),seuil_neb)*rain_int_min))) ELSE znebprecipclr(i)=0.0 ENDIF IF ((zrflcld(i) + ziflcld(i)) .GT. 0.) THEN znebprecipcld(i) = min(znebprecipcld(i), max( & zrflcld(i)/ (MAX(znebprecipcld(i),seuil_neb)*rain_int_min), & ziflcld(i)/ (MAX(znebprecipcld(i),seuil_neb)*rain_int_min))) ELSE znebprecipcld(i)=0.0 ENDIF ENDDO ENDIF ! we recalculate zradoice to account for contributions from the precipitation flux ! if ok_radocond_snow is true IF ( ok_radocond_snow ) THEN IF ( .NOT. iftop ) THEN ! for the solid phase (crystals + snowflakes) ! we recalculate zradoice by summing ! the ice content calculated in the mesh ! + the contribution of the non-evaporated snowfall ! from the overlying layer DO i=1,klon IF (ziflprev(i) .NE. 0.0) THEN zradoice(i)=zoliqi(i)+zqpreci(i)+ziflprev(i)/zrho_up(i)/zvelo_up(i) ELSE zradoice(i)=zoliqi(i)+zqpreci(i) ENDIF zradocond(i)=zradoliq(i)+zradoice(i) ENDDO ENDIF ! keep in memory air density and ice fallspeed velocity zrho_up(:) = zrho(:) zvelo_up(:) = zvelo(:) ENDIF END SUBROUTINE histprecip_postcld !---------------------------------------------------------------- ! Computes the processes-oriented precipitation formulations for ! evaporation and sublimation ! SUBROUTINE poprecip_precld( & klon, dtime, iftop, paprsdn, paprsup, pplay, temp, tempupnew, qvap, & qprecip, precipfracclr, precipfraccld, qvapclrup, qtotupnew, & cldfra, rvc_seri, qliq, qice, & rain, rainclr, raincld, snow, snowclr, snowcld, dqreva, dqssub & ) USE lmdz_lscp_ini, ONLY : prt_level, lunout USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub, expo_eva, expo_sub, thresh_precip_frac USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG USE lmdz_lscp_ini, ONLY : ok_corr_vap_evasub, ok_ice_supersat, ok_unadjusted_clouds USE lmdz_lscp_ini, ONLY : eps, temp_nowater, ok_growth_precip_deposition 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) :: iftop !--if top of the column 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) :: pplay !--pressure in the middle of the layer [Pa] REAL, INTENT(INOUT), DIMENSION(klon) :: temp !--current temperature [K] REAL, INTENT(IN), DIMENSION(klon) :: tempupnew !--updated temperature of the overlying layer [K] REAL, INTENT(INOUT), DIMENSION(klon) :: qvap !--current water vapor specific humidity (includes evaporated qi and ql) [kg/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: qprecip !--specific humidity in the precipitation falling from the upper layer [kg/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: precipfracclr !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-] REAL, INTENT(INOUT), DIMENSION(klon) :: precipfraccld !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-] REAL, INTENT(IN), DIMENSION(klon) :: qvapclrup !--clear-sky specific humidity IN THE LAYER ABOVE [kg/kg] REAL, INTENT(IN), DIMENSION(klon) :: qtotupnew !--total specific humidity IN THE LAYER ABOVE [kg/kg] REAL, INTENT(IN), DIMENSION(klon) :: cldfra !--cloud fraction at the beginning of lscp - used only if the cloud properties are advected [-] REAL, INTENT(IN), DIMENSION(klon) :: rvc_seri !--cloud water vapor at the beginning of lscp (ratio wrt total water) - used only if the cloud properties are advected [kg/kg] REAL, INTENT(IN), DIMENSION(klon) :: qliq !--liquid water content at the beginning of lscp - used only if the cloud properties are advected [kg/kg] REAL, INTENT(IN), DIMENSION(klon) :: qice !--ice water content at the beginning of lscp - used only if the cloud properties are advected [kg/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: rain !--flux of rain gridbox-mean coming from the layer above [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: rainclr !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: raincld !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: snow !--flux of snow gridbox-mean coming from the layer above [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: snowclr !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: snowcld !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2] REAL, INTENT(OUT), DIMENSION(klon) :: dqreva !--rain tendency due to evaporation [kg/kg/s] REAL, INTENT(OUT), DIMENSION(klon) :: dqssub !--snow tendency due to sublimation [kg/kg/s] !--Integer for interating over klon INTEGER :: i !--dhum_to_dflux: coef to convert a specific quantity variation to a flux variation REAL, DIMENSION(klon) :: dhum_to_dflux !--Air density [kg/m3] and layer thickness [m] REAL, DIMENSION(klon) :: rho, dz !--Temporary precip fractions and fluxes for the evaporation REAL, DIMENSION(klon) :: precipfracclr_tmp, precipfraccld_tmp REAL, DIMENSION(klon) :: rainclr_tmp, raincld_tmp, snowclr_tmp, snowcld_tmp !--Saturation values REAL, DIMENSION(klon) :: qzero, qsat, dqsat, qsatl, dqsatl, qsati, dqsati !--Vapor in the clear sky and cloud REAL :: qvapclr, qvapcld !--Fluxes tendencies because of evaporation and sublimation REAL :: dprecip_evasub_max, dprecip_evasub_tot REAL :: drainclreva, draincldeva, dsnowclrsub, dsnowcldsub !--Specific humidity tendencies because of evaporation and sublimation REAL :: dqrevap, dqssubl !--Specific heat constant REAL :: cpair, cpw !--Initialisation qzero(:) = 0. dqreva(:) = 0. dqssub(:) = 0. !-- dhum_to_dflux = rho * dz/dt = 1 / g * dP/dt dhum_to_dflux(:) = ( paprsdn(:) - paprsup(:) ) / RG / dtime rho(:) = pplay(:) / temp(:) / RD dz(:) = ( paprsdn(:) - paprsup(:) ) / RG / rho(:) !--Calculation of saturation specific humidity !--depending on temperature: CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,0,.false.,qsat(:),dqsat(:)) !--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(:)) !--First step consists in "thermalizing" the layer: !--as the flux of precip from layer above "advects" some heat (as the precip is at the temperature !--of the overlying layer) we recalculate a mean temperature that both the air and the precip in the !--layer have. IF (iftop) THEN DO i = 1, klon qprecip(i) = 0. ENDDO ELSE DO i = 1, klon !--No condensed water so cp=cp(vapor+dry air) !-- RVTMP2=rcpv/rcpd-1 cpair = RCPD * ( 1. + RVTMP2 * qvap(i) ) cpw = RCPD * RVTMP2 !--qprecip has to be thermalized with !--layer's air so that precipitation at the ground has the !--same temperature as the lowermost layer !--we convert the flux into a specific quantity qprecip qprecip(i) = ( rain(i) + snow(i) ) / dhum_to_dflux(i) !-- t(i,k+1) + d_t(i,k+1): new temperature of the overlying layer temp(i) = ( tempupnew(i) * qprecip(i) * cpw + cpair * temp(i) ) & / ( cpair + qprecip(i) * cpw ) ENDDO ENDIF !--Initialise the precipitation fractions and fluxes DO i = 1, klon precipfracclr_tmp(i) = precipfracclr(i) precipfraccld_tmp(i) = precipfraccld(i) rainclr_tmp(i) = rainclr(i) snowclr_tmp(i) = snowclr(i) raincld_tmp(i) = raincld(i) snowcld_tmp(i) = snowcld(i) ENDDO !--If we relax the LTP assumption that the cloud is the same than in the !--gridbox above, we can change the tmp quantities IF ( ok_ice_supersat ) THEN !--Update the precipitation fraction using advected cloud fraction CALL poprecip_fracupdate( & klon, cldfra, precipfracclr_tmp, precipfraccld_tmp, & rainclr_tmp, raincld_tmp, snowclr_tmp, snowcld_tmp) ! here we could put an ELSE statement and do one iteration of the condensation scheme ! (but it would only diagnose cloud from lognormal - link with thermals?) DO i = 1, klon !--We cannot have a total fraction greater than the one in the gridbox above !--because new cloud formation has not yet occured. !--If this happens, we reduce the precip frac in the cloud, because it can !--only mean that the cloud fraction at this level is greater than the total !--precipitation fraction of the level above, and the precip frac in clear sky !--is zero. !--NB. this only works because we assume a max-random cloud and precip overlap precipfraccld_tmp(i) = MIN( precipfraccld_tmp(i), precipfracclr(i) + precipfraccld(i) ) ENDDO ENDIF !--At this stage, we guarantee that !-- precipfracclr + precipfraccld = precipfracclr_tmp + precipfraccld_tmp DO i = 1, klon dqrevap = 0. dqssubl = 0. !--If there is precipitation from the layer above IF ( ( rain(i) + snow(i) ) .GT. 0. ) THEN !--Init drainclreva = 0. draincldeva = 0. dsnowclrsub = 0. dsnowcldsub = 0. !--Init the properties of the air where reevaporation occurs IF ( ok_ice_supersat ) THEN !--We can diagnose the water vapor in the cloud using the advected !--water vapor in the cloud and cloud fraction IF ( ok_unadjusted_clouds .AND. ( cldfra(i) .GT. eps ) ) THEN qvapcld = rvc_seri(i) * qvap(i) / cldfra(i) ENDIF !--We can diagnose completely the water vapor in clear sky, because all !--the needed variables (ice, liq, vapor in cld, cloud fraction) are !--advected !--Note that qvap(i) is the total water in the gridbox IF ( ( 1. - cldfra(i) ) .GT. eps ) THEN qvapclr = ( qvap(i) - qice(i) - qliq(i) - rvc_seri(i) * qvap(i) ) / ( 1. - cldfra(i) ) ENDIF ELSEIF ( ok_corr_vap_evasub ) THEN !--Corrected version from Ludo - we use the same water ratio between !--the clear and the cloudy sky as in the layer above. This !--extends the assumption that the cloud fraction is the same !--as the layer above. This is assumed only for the evap / subl !--process !--Note that qvap(i) is the total water in the gridbox, and !--precipfraccld(i) is the cloud fraction in the layer above IF ( ( 1. - precipfraccld(i) ) .GT. eps ) THEN qvapclr = qvapclrup(i) / qtotupnew(i) * qvap(i) / ( 1. - precipfraccld(i) ) ELSE qvapclr = 0. ENDIF IF ( precipfraccld(i) .GT. eps ) THEN qvapcld = MAX(qtotupnew(i)-qvapclrup(i) , 0.) / qtotupnew(i) * qvap(i) / precipfraccld(i) ELSE qvapcld = 0. ENDIF ELSE !--Legacy version from Ludo - we use the total specific humidity !--for the evap / subl process qvapclr = qvap(i) qvapcld = qvap(i) ENDIF !--------------------------- !-- EVAP / SUBL IN CLEAR SKY !--------------------------- IF ( precipfracclr_tmp(i) .GT. eps ) THEN !--Evaporation of liquid precipitation coming from above !--in the clear sky only !--dprecip/dz = -beta*(1-qvap/qsat)*(precip**expo_eva) !--formula from Sundqvist 1988, Klemp & Wilhemson 1978 !--Exact explicit formulation (rainclr is resolved exactly, qvap explicitly) !--which does not need a barrier on rainclr, because included in the formula drainclreva = precipfracclr_tmp(i) * MAX(0., & - coef_eva * ( 1. - expo_eva ) * (1. - qvapclr / qsatl(i)) * dz(i) & + ( rainclr_tmp(i) / precipfracclr_tmp(i) )**( 1. - expo_eva ) & )**( 1. / ( 1. - expo_eva ) ) - rainclr_tmp(i) !--Evaporation is limited by 0 !--NB. with ok_ice_supersat activated, this barrier should be useless drainclreva = MIN(0., drainclreva) ! we set it to 0 as not sufficiently tested drainclreva = 0. !--Sublimation of the solid precipitation coming from above !--(same formula as for liquid precip) !--Exact explicit formulation (snowclr is resolved exactly, qvap explicitly) !--which does not need a barrier on snowclr, because included in the formula dsnowclrsub = precipfracclr_tmp(i) * MAX(0., & - coef_sub * ( 1. - expo_sub ) * (1. - qvapclr / qsati(i)) * dz(i) & + ( snowclr_tmp(i) / precipfracclr_tmp(i) )**( 1. - expo_sub ) & )**( 1. / ( 1. - expo_sub ) ) - snowclr_tmp(i) !--If ice supersaturation is activated, we allow vapor to condense on falling snow !--i.e., having a positive dsnowclrsub IF ( .NOT. ok_ice_supersat ) THEN !--Sublimation is limited by 0 dsnowclrsub = MIN(0., dsnowclrsub) ENDIF !--Evaporation limit: we ensure that the layer's fraction below !--the clear sky does not reach saturation. In this case, we !--redistribute the maximum flux dprecip_evasub_max conserving the ratio liquid/ice !--Max evaporation is computed not to saturate the clear sky precip fraction !--(i.e., the fraction where evaporation occurs) !--It is expressed as a max flux dprecip_evasub_max IF ( .NOT. ok_ice_supersat ) THEN dprecip_evasub_max = MIN(0., ( qvapclr - qsat(i) ) * precipfracclr_tmp(i)) & * dhum_to_dflux(i) ELSE dprecip_evasub_max = ( qvapclr - qsat(i) ) * precipfracclr_tmp(i) & * dhum_to_dflux(i) ENDIF dprecip_evasub_tot = drainclreva + dsnowclrsub !--Barriers !--If activates if the total is LOWER than the max because !--everything is negative IF ( (( dprecip_evasub_max .LT. 0. ) .AND. & ( dprecip_evasub_tot .LT. dprecip_evasub_max )) .OR. & (( dprecip_evasub_max .GT. 0. ) .AND. & ( dprecip_evasub_tot .GT. dprecip_evasub_max )) ) THEN drainclreva = dprecip_evasub_max * drainclreva / dprecip_evasub_tot dsnowclrsub = dprecip_evasub_max * dsnowclrsub / dprecip_evasub_tot ENDIF ELSE !--All the precipitation is sublimated if the fraction is zero drainclreva = - rainclr_tmp(i) dsnowclrsub = - snowclr_tmp(i) ENDIF ! precipfracclr_tmp .GT. eps !--------------------------- !-- EVAP / SUBL IN THE CLOUD !--------------------------- IF ( ok_growth_precip_deposition .AND. ( temp(i) .LE. RTT ) .AND. ( precipfraccld_tmp(i) .GT. eps ) ) THEN !--Evaporation of liquid precipitation coming from above !--in the cloud only !--dprecip/dz = -beta*(1-qvap/qsat)*(precip**expo_eva) !--formula from Sundqvist 1988, Klemp & Wilhemson 1978 !--Exact explicit formulation (raincld is resolved exactly, qvap explicitly) !--which does not need a barrier on raincld, because included in the formula draincldeva = precipfraccld_tmp(i) * MAX(0., & - coef_eva * ( 1. - expo_eva ) * (1. - qvapcld / qsatl(i)) * dz(i) & + ( raincld_tmp(i) / precipfraccld_tmp(i) )**( 1. - expo_eva ) & )**( 1. / ( 1. - expo_eva ) ) - raincld_tmp(i) !--Evaporation is limited by 0 !--NB. with ok_ice_supersat activated, this barrier should be useless draincldeva = MIN(0., draincldeva) !--Sublimation of the solid precipitation coming from above !--(same formula as for liquid precip) !--Exact explicit formulation (snowcld is resolved exactly, qvap explicitly) !--which does not need a barrier on snowcld, because included in the formula dsnowcldsub = precipfraccld_tmp(i) * MAX(0., & - coef_sub * ( 1. - expo_sub ) * (1. - qvapcld / qsati(i)) * dz(i) & + ( snowcld_tmp(i) / precipfraccld_tmp(i) )**( 1. - expo_sub ) & )**( 1. / ( 1. - expo_sub ) ) - snowcld_tmp(i) !--There is no barrier because we want deposition to occur if it is possible !--Evaporation limit: we ensure that the layer's fraction below !--the clear sky does not reach saturation. In this case, we !--redistribute the maximum flux dprecip_evasub_max conserving the ratio liquid/ice !--Max evaporation is computed not to saturate the clear sky precip fraction !--(i.e., the fraction where evaporation occurs) !--It is expressed as a max flux dprecip_evasub_max dprecip_evasub_max = ( qvapcld - qsat(i) ) * precipfraccld_tmp(i) * dhum_to_dflux(i) dprecip_evasub_tot = draincldeva + dsnowcldsub !--Barriers !--If activates if the total is LOWER than the max because !--everything is negative IF ( (( dprecip_evasub_max .LT. 0. ) .AND. & ( dprecip_evasub_tot .LT. dprecip_evasub_max )) .OR. & (( dprecip_evasub_max .GT. 0. ) .AND. & ( dprecip_evasub_tot .GT. dprecip_evasub_max )) ) THEN draincldeva = dprecip_evasub_max * draincldeva / dprecip_evasub_tot dsnowcldsub = dprecip_evasub_max * dsnowcldsub / dprecip_evasub_tot ENDIF ENDIF ! ok_unadjusted_clouds !--New solid and liquid precipitation fluxes after evap and sublimation dqrevap = ( drainclreva + draincldeva ) / dhum_to_dflux(i) dqssubl = ( dsnowclrsub + dsnowcldsub ) / dhum_to_dflux(i) !--Vapor is updated after evaporation/sublimation (it is increased) qvap(i) = qvap(i) - dqrevap - dqssubl !--qprecip is the total condensed water in the precip flux (it is decreased) qprecip(i) = qprecip(i) + dqrevap + dqssubl !--Air and precip temperature (i.e., gridbox temperature) !--is updated due to latent heat cooling temp(i) = temp(i) & + dqrevap * RLVTT / RCPD & / ( 1. + RVTMP2 * ( qvap(i) + qprecip(i) ) ) & + dqssubl * RLSTT / RCPD & / ( 1. + RVTMP2 * ( qvap(i) + qprecip(i) ) ) !--Add tendencies !--The MAX is needed because in some cases, the flux can be slightly !--negative (numerical precision) rainclr_tmp(i) = MAX(0., rainclr_tmp(i) + drainclreva) snowclr_tmp(i) = MAX(0., snowclr_tmp(i) + dsnowclrsub) raincld_tmp(i) = MAX(0., raincld_tmp(i) + draincldeva) snowcld_tmp(i) = MAX(0., snowcld_tmp(i) + dsnowcldsub) !--We reattribute fluxes according to initial fractions, using proportionnality !--Note that trough this process, we conserve total precip fluxes !-- rainclr_tmp + raincld_tmp = rainclr + raincld (same for snow) !--If the cloud has increased, a part of the precip in clear sky falls in clear sky, !--the rest falls in the cloud IF ( ( precipfraccld(i) .GT. eps ) .AND. ( precipfraccld_tmp(i) .GT. precipfraccld(i) ) ) THEN !--All the precip from cloud falls in the cloud raincld(i) = MAX(0., raincld_tmp(i) * precipfraccld(i) / precipfraccld_tmp(i)) rainclr(i) = MAX(0., rainclr_tmp(i) + raincld_tmp(i) - raincld(i)) snowcld(i) = MAX(0., snowcld_tmp(i) * precipfraccld(i) / precipfraccld_tmp(i)) snowclr(i) = MAX(0., snowclr_tmp(i) + snowcld_tmp(i) - snowcld(i)) !--If the cloud has narrowed, a part of the precip in cloud falls in clear sky, !--the rest falls in the cloud ELSEIF ( ( precipfracclr(i) .GT. eps ) .AND. ( precipfracclr_tmp(i) .GT. precipfracclr(i) ) ) THEN !--All the precip from clear sky falls in clear sky rainclr(i) = MAX(0., rainclr_tmp(i) * precipfracclr(i) / precipfracclr_tmp(i)) raincld(i) = MAX(0., raincld_tmp(i) + rainclr_tmp(i) - rainclr(i)) snowclr(i) = MAX(0., snowclr_tmp(i) * precipfracclr(i) / precipfracclr_tmp(i)) snowcld(i) = MAX(0., snowcld_tmp(i) + snowclr_tmp(i) - snowclr(i)) ELSE !--If the cloud stays the same or if there is no cloud above and !--in the current layer, there is no reattribution rainclr(i) = rainclr_tmp(i) raincld(i) = raincld_tmp(i) snowclr(i) = snowclr_tmp(i) snowcld(i) = snowcld_tmp(i) ENDIF !--If there is no more precip fluxes, the precipitation fraction is set to 0 IF ( ( rainclr(i) + snowclr(i) ) .LE. 0. ) precipfracclr(i) = 0. IF ( ( raincld(i) + snowcld(i) ) .LE. 0. ) precipfraccld(i) = 0. !--Calculation of the total fluxes rain(i) = rainclr(i) + raincld(i) snow(i) = snowclr(i) + snowcld(i) ELSEIF ( ( rain(i) + snow(i) ) .LE. 0. ) THEN !--If no precip, we reinitialize the cloud fraction used for the precip to 0 precipfraccld(i) = 0. precipfracclr(i) = 0. ENDIF ! ( ( rain(i) + snow(i) ) .GT. 0. ) !--Diagnostic tendencies dqssub(i) = dqssubl / dtime dqreva(i) = dqrevap / dtime ENDDO ! loop on klon END SUBROUTINE poprecip_precld !---------------------------------------------------------------- ! Computes the precipitation fraction update ! The goal of this routine is to reattribute precipitation fractions ! and fluxes to clear or cloudy air, depending on the variation of ! the cloud fraction on the vertical dimension. We assume a ! maximum-random overlap of the cloud cover (see Jakob and Klein, 2000, ! and LTP thesis, 2021) ! NB. in fact, we assume a maximum-random overlap of the total precip. frac ! SUBROUTINE poprecip_fracupdate( & klon, cldfra, precipfracclr, precipfraccld, & rainclr, raincld, snowclr, snowcld) USE lmdz_lscp_ini, ONLY : eps IMPLICIT NONE INTEGER, INTENT(IN) :: klon !--number of horizontal grid points [-] REAL, INTENT(IN), DIMENSION(klon) :: cldfra !--cloud fraction [-] REAL, INTENT(INOUT), DIMENSION(klon) :: precipfracclr !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-] REAL, INTENT(INOUT), DIMENSION(klon) :: precipfraccld !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-] !--NB. at the end of the routine, becomes the fraction of precip !--in the current layer REAL, INTENT(INOUT), DIMENSION(klon) :: rainclr !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: raincld !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: snowclr !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: snowcld !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2] !--Local variables INTEGER :: i REAL :: dcldfra REAL :: precipfractot REAL :: dprecipfracclr, dprecipfraccld REAL :: drainclr, dsnowclr REAL :: draincld, dsnowcld DO i = 1, klon !--Initialisation precipfractot = precipfracclr(i) + precipfraccld(i) !--Instead of using the cloud cover which was use in LTP thesis, we use the !--total precip. fraction to compute the maximum-random overlap. This is !--because all the information of the cloud cover is embedded into !--precipfractot, and this allows for taking into account the potential !--reduction of the precipitation fraction because either the flux is too !--small (see barrier at the end of poprecip_postcld) or the flux is completely !--evaporated (see barrier at the end of poprecip_precld) !--NB. precipfraccld(i) is here the cloud fraction of the layer above !precipfractot = 1. - ( 1. - precipfractot ) * & ! ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) & ! / ( 1. - MIN( precipfraccld(i), 1. - eps ) ) IF ( precipfraccld(i) .GT. ( 1. - eps ) ) THEN precipfractot = 1. ELSE precipfractot = 1. - ( 1. - precipfractot ) * & ( 1. - MAX( cldfra(i), precipfraccld(i) ) ) & / ( 1. - precipfraccld(i) ) ENDIF !--precipfraccld(i) is here the cloud fraction of the layer above dcldfra = cldfra(i) - precipfraccld(i) !--Tendency of the clear-sky precipitation fraction. We add a MAX on the !--calculation of the current CS precip. frac. !dprecipfracclr = MAX( 0., ( precipfractot - cldfra(i) ) ) - precipfracclr(i) !--We remove it, because precipfractot is guaranteed to be > cldfra (the MAX is activated !--if precipfractot < cldfra) dprecipfracclr = ( precipfractot - cldfra(i) ) - precipfracclr(i) !--Tendency of the cloudy precipitation fraction. We add a MAX on the !--calculation of the current CS precip. frac. !dprecipfraccld = MAX( dcldfra , - precipfraccld(i) ) !--We remove it, because cldfra is guaranteed to be > 0 (the MAX is activated !--if cldfra < 0) dprecipfraccld = dcldfra !--If the cloud extends IF ( dprecipfraccld .GT. 0. ) THEN !--If there is no CS precip, nothing happens. !--If there is, we reattribute some of the CS precip flux !--to the cloud precip flux, proportionnally to the !--decrease of the CS precip fraction IF ( precipfracclr(i) .LT. eps ) THEN drainclr = 0. dsnowclr = 0. ELSE drainclr = dprecipfracclr / precipfracclr(i) * rainclr(i) dsnowclr = dprecipfracclr / precipfracclr(i) * snowclr(i) ENDIF !--If the cloud narrows ELSEIF ( dprecipfraccld .LT. 0. ) THEN !--We reattribute some of the cloudy precip flux !--to the CS precip flux, proportionnally to the !--decrease of the cloud precip fraction IF ( precipfraccld(i) .LT. eps ) THEN draincld = 0. dsnowcld = 0. ELSE draincld = dprecipfraccld / precipfraccld(i) * raincld(i) dsnowcld = dprecipfraccld / precipfraccld(i) * snowcld(i) ENDIF drainclr = - draincld dsnowclr = - dsnowcld !--If the cloud stays the same or if there is no cloud above and !--in the current layer, nothing happens ELSE drainclr = 0. dsnowclr = 0. ENDIF !--We add the tendencies !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision) precipfraccld(i) = precipfraccld(i) + dprecipfraccld precipfracclr(i) = precipfracclr(i) + dprecipfracclr rainclr(i) = MAX(0., rainclr(i) + drainclr) snowclr(i) = MAX(0., snowclr(i) + dsnowclr) raincld(i) = MAX(0., raincld(i) - drainclr) snowcld(i) = MAX(0., snowcld(i) - dsnowclr) ENDDO END SUBROUTINE poprecip_fracupdate !---------------------------------------------------------------- ! Computes the processes-oriented precipitation formulations for ! - autoconversion (auto) via a deposition process ! - aggregation (agg) ! - riming (rim) ! - collection (col) ! - melting (melt) ! - freezing (freez) ! SUBROUTINE poprecip_postcld( & klon, dtime, paprsdn, paprsup, pplay, ctot_vol, ptconv, & temp, qvap, qliq, qice, icefrac, cldfra, & precipfracclr, precipfraccld, & rain, rainclr, raincld, snow, snowclr, snowcld, & qraindiag, qsnowdiag, dqrauto, dqrcol, dqrmelt, dqrfreez, & dqsauto, dqsagg, dqsrim, dqsmelt, dqsfreez) USE lmdz_lscp_ini, ONLY : prt_level, lunout USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RPI USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf USE lmdz_lscp_ini, ONLY : cld_lc_con, cld_tau_con, cld_expo_con, seuil_neb, & cld_lc_lsc, cld_tau_lsc, cld_expo_lsc, & thresh_precip_frac, gamma_col, gamma_agg, gamma_rim, & rho_rain, r_rain, r_snow, rho_ice, & expo_tau_auto_snow, & tau_auto_snow_min, tau_auto_snow_max, & thresh_precip_frac, eps, rain_int_min, & gamma_melt, alpha_freez, beta_freez, temp_nowater, & iflag_cloudth_vert, iflag_rain_incloud_vol, & cld_lc_lsc_snow, cld_lc_con_snow, gamma_freez, & rain_fallspeed_clr, rain_fallspeed_cld, & snow_fallspeed_clr, snow_fallspeed_cld IMPLICIT NONE INTEGER, INTENT(IN) :: klon !--number of horizontal grid points [-] REAL, INTENT(IN) :: dtime !--time step [s] 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) :: pplay !--pressure in the middle of the layer [Pa] REAL, INTENT(IN), DIMENSION(klon) :: ctot_vol !--volumic cloud fraction [-] LOGICAL, INTENT(IN), DIMENSION(klon) :: ptconv !--true if we are in a convective point REAL, INTENT(INOUT), DIMENSION(klon) :: temp !--current temperature [K] REAL, INTENT(INOUT), DIMENSION(klon) :: qvap !--current water vapor specific humidity [kg/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: qliq !--current liquid water specific humidity [kg/kg] REAL, INTENT(INOUT), DIMENSION(klon) :: qice !--current ice water specific humidity [kg/kg] REAL, INTENT(IN), DIMENSION(klon) :: icefrac !--ice fraction [-] REAL, INTENT(IN), DIMENSION(klon) :: cldfra !--cloud fraction [-] REAL, INTENT(INOUT), DIMENSION(klon) :: precipfracclr !--fraction of precipitation in the clear sky IN THE LAYER ABOVE [-] REAL, INTENT(INOUT), DIMENSION(klon) :: precipfraccld !--fraction of precipitation in the cloudy air IN THE LAYER ABOVE [-] !--NB. at the end of the routine, becomes the fraction of precip !--in the current layer REAL, INTENT(INOUT), DIMENSION(klon) :: rain !--flux of rain gridbox-mean coming from the layer above [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: rainclr !--flux of rain gridbox-mean in clear sky coming from the layer above [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: raincld !--flux of rain gridbox-mean in cloudy air coming from the layer above [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: snow !--flux of snow gridbox-mean coming from the layer above [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: snowclr !--flux of snow gridbox-mean in clear sky coming from the layer above [kg/s/m2] REAL, INTENT(INOUT), DIMENSION(klon) :: snowcld !--flux of snow gridbox-mean in cloudy air coming from the layer above [kg/s/m2] REAL, INTENT(OUT), DIMENSION(klon) :: qraindiag !--DIAGNOSTIC specific rain content [kg/kg] REAL, INTENT(OUT), DIMENSION(klon) :: qsnowdiag !--DIAGNOSTIC specific snow content [kg/kg] REAL, INTENT(OUT), DIMENSION(klon) :: dqrcol !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s] REAL, INTENT(OUT), DIMENSION(klon) :: dqsagg !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s] REAL, INTENT(OUT), DIMENSION(klon) :: dqrauto !--rain tendency due to autoconversion of cloud liquid [kg/kg/s] REAL, INTENT(OUT), DIMENSION(klon) :: dqsauto !--snow tendency due to autoconversion of cloud ice [kg/kg/s] REAL, INTENT(OUT), DIMENSION(klon) :: dqsrim !--snow tendency due to riming [kg/kg/s] REAL, INTENT(OUT), DIMENSION(klon) :: dqsmelt !--snow tendency due to melting [kg/kg/s] REAL, INTENT(OUT), DIMENSION(klon) :: dqrmelt !--rain tendency due to melting [kg/kg/s] REAL, INTENT(OUT), DIMENSION(klon) :: dqsfreez !--snow tendency due to freezing [kg/kg/s] REAL, INTENT(OUT), DIMENSION(klon) :: dqrfreez !--rain tendency due to freezing [kg/kg/s] !--Local variables INTEGER :: i !--dhum_to_dflux: coef to convert a specific quantity variation to a flux variation REAL, DIMENSION(klon) :: dhum_to_dflux REAL, DIMENSION(klon) :: qtot !--includes vap, liq, ice and precip REAL :: min_precip !--minimum precip flux below which precip fraction decreases !--Collection, aggregation and riming REAL :: eff_cldfra REAL :: coef_col, coef_agg, coef_rim, coef_tmp, qrain_tmp REAL :: Eff_rain_liq, Eff_snow_ice, Eff_snow_liq REAL :: rho_snow REAL :: dqlcol !--loss of liquid cloud content due to collection by rain [kg/kg/s] REAL :: dqiagg !--loss of ice cloud content due to collection by aggregation [kg/kg/s] REAL :: dqlrim !--loss of liquid cloud content due to riming on snow [kg/kg/s] !--Autoconversion REAL :: qthresh_auto_rain, tau_auto_rain, expo_auto_rain REAL :: qthresh_auto_snow, tau_auto_snow, expo_auto_snow REAL :: dqlauto !--loss of liquid cloud content due to autoconversion to rain [kg/kg/s] REAL :: dqiauto !--loss of ice cloud content due to autoconversion to snow [kg/kg/s] !--Melting REAL :: dqsmelt_max, air_thermal_conduct REAL :: nb_snowflake_clr, nb_snowflake_cld REAL :: capa_snowflake, temp_wetbulb REAL :: rho, r_ice REAL :: dqsclrmelt, dqscldmelt, dqstotmelt REAL, DIMENSION(klon) :: qzero, qsat, dqsat !--Freezing REAL :: dqrfreez_max REAL :: tau_freez REAL :: dqrclrfreez, dqrcldfreez, dqrtotfreez, dqrtotfreez_step1, dqrtotfreez_step2 REAL :: coef_freez REAL :: dqifreez !--loss of ice cloud content due to collection of ice from rain [kg/kg/s] REAL :: Eff_rain_ice !--Initialisation of variables qzero(:) = 0. dqrcol(:) = 0. dqsagg(:) = 0. dqsauto(:) = 0. dqrauto(:) = 0. dqsrim(:) = 0. dqrmelt(:) = 0. dqsmelt(:) = 0. dqrfreez(:) = 0. dqsfreez(:) = 0. !--Update the precipitation fraction following cloud formation CALL poprecip_fracupdate( & klon, cldfra, precipfracclr, precipfraccld, & rainclr, raincld, snowclr, snowcld) DO i = 1, klon !--Variables initialisation dqlcol = 0. dqiagg = 0. dqiauto = 0. dqlauto = 0. dqlrim = 0. !-- dhum_to_dflux = rho * dz/dt = 1 / g * dP/dt dhum_to_dflux(i) = ( paprsdn(i) - paprsup(i) ) / RG / dtime qtot(i) = qvap(i) + qliq(i) + qice(i) & + ( raincld(i) + rainclr(i) + snowcld(i) + snowclr(i) ) / dhum_to_dflux(i) !--If vertical heterogeneity is taken into account, we use !--the "true" volume fraction instead of a modified !--surface fraction (which is larger and artificially !--reduces the in-cloud water). IF ( ( iflag_cloudth_vert .GE. 3 ) .AND. ( iflag_rain_incloud_vol .EQ. 1 ) ) THEN eff_cldfra = ctot_vol(i) ELSE eff_cldfra = cldfra(i) ENDIF !--Start precipitation growth processes !--If the cloud is big enough, the precipitation processes activate ! TODO met on seuil_neb ici ? IF ( cldfra(i) .GE. seuil_neb ) THEN !--------------------------------------------------------- !-- COLLECTION AND AGGREGATION !--------------------------------------------------------- !--Collection: processus through which rain collects small liquid droplets !--in suspension, and add it to the rain flux !--Aggregation: same for snow (precip flux) and ice crystals (in suspension) !--Those processes are treated before autoconversion because we do not !--want to collect/aggregate the newly formed fluxes, which already !--"saw" the cloud as they come from it !--The formulas come from Muench and Lohmann 2020 !--gamma_col: tuning coefficient [-] !--rho_rain: volumic mass of rain [kg/m3] !--r_rain: size of the rain droplets [m] !--Eff_rain_liq: efficiency of the collection process [-] (between 0 and 1) !--dqlcol is a gridbox-mean quantity, as is qliq and raincld. They are !--divided by respectively eff_cldfra, eff_cldfra and precipfraccld to !--get in-cloud mean quantities. The two divisions by eff_cldfra are !--then simplified. !--The collection efficiency is perfect. Eff_rain_liq = 1. coef_col = gamma_col * 3. / 4. / rho_rain / r_rain * Eff_rain_liq IF ( raincld(i) .GT. 0. ) THEN !--Exact explicit version, which does not need a barrier because of !--the exponential decrease dqlcol = qliq(i) * ( EXP( - dtime * coef_col * raincld(i) / precipfraccld(i) ) - 1. ) !--Add tendencies qliq(i) = qliq(i) + dqlcol raincld(i) = raincld(i) - dqlcol * dhum_to_dflux(i) !--Diagnostic tendencies dqrcol(i) = - dqlcol / dtime ENDIF !--Same as for aggregation !--Eff_snow_liq formula: !--it s a product of a collection efficiency and a sticking efficiency ! Milbrandt and Yau formula that gives very low values: ! Eff_snow_ice = 0.05 * EXP( 0.1 * ( temp(i) - RTT ) ) ! Lin 1983's formula Eff_snow_ice = EXP( 0.025 * MIN( ( temp(i) - RTT ), 0.) ) !--rho_snow formula follows Brandes et al. 2007 (JAMC) rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922) coef_agg = gamma_agg * 3. / 4. / rho_snow / r_snow * Eff_snow_ice IF ( snowcld(i) .GT. 0. ) THEN !--Exact explicit version, which does not need a barrier because of !--the exponential decrease dqiagg = qice(i) * ( EXP( - dtime * coef_agg * snowcld(i) / precipfraccld(i) ) - 1. ) !--Add tendencies qice(i) = qice(i) + dqiagg snowcld(i) = snowcld(i) - dqiagg * dhum_to_dflux(i) !--Diagnostic tendencies dqsagg(i) = - dqiagg / dtime ENDIF !--------------------------------------------------------- !-- AUTOCONVERSION !--------------------------------------------------------- !--Autoconversion converts liquid droplets/ice crystals into !--rain drops/snowflakes. It relies on the formulations by !--Sundqvist 1978. !--If we are in a convective point, we have different parameters !--for the autoconversion IF ( ptconv(i) ) THEN qthresh_auto_rain = cld_lc_con qthresh_auto_snow = cld_lc_con_snow tau_auto_rain = cld_tau_con !--tau for snow depends on the ice fraction in mixed-phase clouds tau_auto_snow = tau_auto_snow_max & + ( tau_auto_snow_min - tau_auto_snow_max ) * ((1. - icefrac(i))**expo_tau_auto_snow) expo_auto_rain = cld_expo_con expo_auto_snow = cld_expo_con ELSE qthresh_auto_rain = cld_lc_lsc qthresh_auto_snow = cld_lc_lsc_snow tau_auto_rain = cld_tau_lsc !--tau for snow depends on the ice fraction in mixed-phase clouds tau_auto_snow = tau_auto_snow_max & + ( tau_auto_snow_min - tau_auto_snow_max ) * ((1. - icefrac(i))**expo_tau_auto_snow) expo_auto_rain = cld_expo_lsc expo_auto_snow = cld_expo_lsc ENDIF ! Liquid water quantity to remove according to (Sundqvist, 1978) ! dqliq/dt = -qliq/tau * ( 1-exp(-(qliqincld/qthresh)**2) ) ! !--And same formula for ice ! !--We first treat the second term (with exponential) in an explicit way !--and then treat the first term (-q/tau) in an exact way dqlauto = - qliq(i) * ( 1. - exp( - dtime / tau_auto_rain * ( 1. - exp( & - ( qliq(i) / eff_cldfra / qthresh_auto_rain ) ** expo_auto_rain ) ) ) ) dqiauto = - qice(i) * ( 1. - exp( - dtime / tau_auto_snow * ( 1. - exp( & - ( qice(i) / eff_cldfra / qthresh_auto_snow ) ** expo_auto_snow ) ) ) ) !--Barriers so that we don't create more rain/snow !--than there is liquid/ice dqlauto = MAX( - qliq(i), dqlauto ) dqiauto = MAX( - qice(i), dqiauto ) !--Add tendencies qliq(i) = qliq(i) + dqlauto qice(i) = qice(i) + dqiauto raincld(i) = raincld(i) - dqlauto * dhum_to_dflux(i) snowcld(i) = snowcld(i) - dqiauto * dhum_to_dflux(i) !--Diagnostic tendencies dqsauto(i) = - dqiauto / dtime dqrauto(i) = - dqlauto / dtime !--------------------------------------------------------- !-- RIMING !--------------------------------------------------------- !--Process which converts liquid droplets in suspension into !--snow because of the collision between !--those and falling snowflakes. !--The formula comes from Muench and Lohmann 2020 !--NB.: this process needs a temperature adjustment !--Eff_snow_liq formula: following Ferrier 1994, !--assuming 1 Eff_snow_liq = 1.0 !--rho_snow formula follows Brandes et al. 2007 (JAMC) rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922) coef_rim = gamma_rim * 3. / 4. / rho_snow / r_snow * Eff_snow_liq IF ( snowcld(i) .GT. 0. ) THEN !--Exact version, which does not need a barrier because of !--the exponential decrease dqlrim = qliq(i) * ( EXP( - dtime * coef_rim * snowcld(i) / precipfraccld(i) ) - 1. ) !--Add tendencies !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision) qliq(i) = qliq(i) + dqlrim snowcld(i) = snowcld(i) - dqlrim * dhum_to_dflux(i) !--Temperature adjustment with the release of latent !--heat because of solid condensation temp(i) = temp(i) - dqlrim * RLMLT / RCPD & / ( 1. + RVTMP2 * qtot(i) ) !--Diagnostic tendencies dqsrim(i) = - dqlrim / dtime ENDIF ENDIF ! cldfra .GE. seuil_neb ENDDO ! loop on klon !--Re-calculation of saturation specific humidity !--because riming changed temperature CALL calc_qsat_ecmwf(klon, temp, qzero, pplay, RTT, 0, .FALSE., qsat, dqsat) DO i = 1, klon !--------------------------------------------------------- !-- MELTING !--------------------------------------------------------- !--Process through which snow melts into rain. !--The formula is homemade. !--NB.: this process needs a temperature adjustment !--dqsmelt_max : maximum snow melting so that temperature !-- stays higher than 273 K [kg/kg] !--capa_snowflake : capacitance of a snowflake, equal to !-- the radius if the snowflake is a sphere [m] !--temp_wetbulb : wet-bulb temperature [K] !--snow_fallspeed : snow fall velocity (in clear/cloudy sky) [m/s] !--air_thermal_conduct : thermal conductivity of the air [J/m/K/s] !--gamma_melt : tuning parameter for melting [-] !--nb_snowflake : number of snowflakes (in clear/cloudy air) [-] IF ( ( snowclr(i) + snowcld(i) ) .GT. 0. ) THEN !--Computed according to !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q dqsmelt_max = MIN(0., ( RTT - temp(i) ) / RLMLT * RCPD & * ( 1. + RVTMP2 * qtot(i) )) !--Initialisation dqsclrmelt = 0. dqscldmelt = 0. !--We assume that the snowflakes are spherical capa_snowflake = r_snow !--Thermal conductivity of the air, empirical formula from Beard and Pruppacher (1971) air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184 !--rho_snow formula follows Brandes et al. 2007 (JAMC) rho_snow = 1.e3 * 0.178 * ( r_snow * 2. * 1000. )**(-0.922) !--In clear air IF ( ( snowclr(i) .GT. 0. ) .AND. ( precipfracclr(i) .GT. 0. ) ) THEN !--Formula for the wet-bulb temperature from ECMWF (IFS) !--The vapor used is the vapor in the clear sky temp_wetbulb = temp(i) & - ( qsat(i) - ( qvap(i) - cldfra(i) * qsat(i) ) / ( 1. - cldfra(i) ) ) & * ( 1329.31 + 0.0074615 * ( pplay(i) - 0.85e5 ) & - 40.637 * ( temp(i) - 275. ) ) !--Calculated according to !-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow nb_snowflake_clr = snowclr(i) / precipfracclr(i) / snow_fallspeed_clr & / ( 4. / 3. * RPI * r_snow**3. * rho_snow ) dqsclrmelt = - nb_snowflake_clr * 4. * RPI * air_thermal_conduct & * capa_snowflake / RLMLT * gamma_melt & * MAX(0., temp_wetbulb - RTT) * dtime !--Barrier to limit the melting flux to the clr snow flux in the mesh dqsclrmelt = MAX( dqsclrmelt , -snowclr(i) / dhum_to_dflux(i)) ENDIF !--In cloudy air IF ( ( snowcld(i) .GT. 0. ) .AND. ( precipfraccld(i) .GT. 0. ) ) THEN !--As the air is saturated, the wet-bulb temperature is equal to the !--temperature temp_wetbulb = temp(i) !--Calculated according to !-- flux = velocity_snowflakes * nb_snowflakes * volume_snowflakes * rho_snow nb_snowflake_cld = snowcld(i) / precipfraccld(i) / snow_fallspeed_cld & / ( 4. / 3. * RPI * r_snow**3. * rho_snow ) dqscldmelt = - nb_snowflake_cld * 4. * RPI * air_thermal_conduct & * capa_snowflake / RLMLT * gamma_melt & * MAX(0., temp_wetbulb - RTT) * dtime !--Barrier to limit the melting flux to the cld snow flux in the mesh dqscldmelt = MAX(dqscldmelt , - snowcld(i) / dhum_to_dflux(i)) ENDIF !--Barrier on temperature. If the total melting flux leads to a !--positive temperature, it is limited to keep temperature above 0 degC. !--It is activated if the total is LOWER than the max !--because everything is negative dqstotmelt = dqsclrmelt + dqscldmelt IF ( dqstotmelt .LT. dqsmelt_max ) THEN !--We redistribute the max melted snow keeping !--the clear/cloud partition of the melted snow dqsclrmelt = dqsmelt_max * dqsclrmelt / dqstotmelt dqscldmelt = dqsmelt_max * dqscldmelt / dqstotmelt dqstotmelt = dqsmelt_max ENDIF !--Add tendencies !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision) rainclr(i) = MAX(0., rainclr(i) - dqsclrmelt * dhum_to_dflux(i)) raincld(i) = MAX(0., raincld(i) - dqscldmelt * dhum_to_dflux(i)) snowclr(i) = MAX(0., snowclr(i) + dqsclrmelt * dhum_to_dflux(i)) snowcld(i) = MAX(0., snowcld(i) + dqscldmelt * dhum_to_dflux(i)) !--Temperature adjustment with the release of latent !--heat because of melting temp(i) = temp(i) + dqstotmelt * RLMLT / RCPD & / ( 1. + RVTMP2 * qtot(i) ) !--Diagnostic tendencies dqrmelt(i) = - dqstotmelt / dtime dqsmelt(i) = dqstotmelt / dtime ENDIF !--------------------------------------------------------- !-- FREEZING !--------------------------------------------------------- !--Process through which rain freezes into snow. !-- We parameterize it as a 2 step process: !--first: freezing following collision with ice crystals !--second: immersion freezing following (inspired by Bigg 1953) !--the latter is parameterized as an exponential decrease of the rain !--water content with a homemade formula !--This is based on a caracteritic time of freezing, which !--exponentially depends on temperature so that it is !--equal to 1 for temp_nowater (see below) and is close to !--0 for RTT (=273.15 K). !--NB.: this process needs a temperature adjustment !--dqrfreez_max : maximum rain freezing so that temperature !-- stays lower than 273 K [kg/kg] !--tau_freez : caracteristic time of freezing [s] !--gamma_freez : tuning parameter [s-1] !--alpha_freez : tuning parameter for the shape of the exponential curve [-] !--temp_nowater : temperature below which no liquid water exists [K] (about -40 degC) IF ( ( rainclr(i) + raincld(i) ) .GT. 0. ) THEN !--1st step: freezing following collision with ice crystals !--Sub-process of freezing which quantifies the collision between !--ice crystals in suspension and falling rain droplets. !--The rain droplets freeze, becoming graupel, and carrying !--the ice crystal (which acted as an ice nucleating particle). !--The formula is adapted from the riming formula. !--it works only in the cloudy part dqifreez = 0. dqrtotfreez_step1 = 0. IF ( ( qice(i) .GT. 0. ) .AND. ( cldfra(i) .GT. 0. ) .AND. & ( raincld(i) .GT. 0. ) .AND. ( precipfraccld(i) .GT. 0. ) ) THEN dqrclrfreez = 0. dqrcldfreez = 0. !--Computed according to !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD & * ( 1. + RVTMP2 * qtot(i) )) !--The collision efficiency is assumed unity Eff_rain_ice = 1. coef_freez = gamma_freez * 3. / 4. / rho_rain / r_rain * Eff_rain_ice !--Exact version, which does not need a barrier because of !--the exponential decrease. dqifreez = qice(i) * ( EXP( - dtime * coef_freez * raincld(i) / precipfraccld(i) ) - 1. ) !--We add the part of rain water that freezes, limited by a temperature barrier !--This quantity is calculated assuming that the number of drop that freeze correspond to the number !--of crystals collected (and assuming uniform distributions of ice crystals and rain drops) !--The ice specific humidity that collide with rain is dqi = dNi 4/3 PI rho_ice r_ice**3 !--The rain that collide with ice is, similarly, dqr = dNr 4/3 PI rho_rain r_rain**3 !--The assumption above corresponds to dNi = dNr, i.e., !-- dqr = dqi * (4/3 PI rho_rain * r_rain**3) / (4/3 PI rho_ice * r_ice**3) !--Dry density [kg/m3] rho = pplay(i) / temp(i) / RD !--r_ice formula from Sun and Rikus (1999) r_ice = 1.e-6 * ( 45.8966 * ( qice(i) / cldfra(i) * rho * 1e3 )**0.2214 & + 0.7957 * ( qice(i) / cldfra(i) * rho * 1e3 )**0.2535 * ( temp(i) - RTT + 190. ) ) / 2. dqrcldfreez = dqifreez * rho_rain * r_rain**3. / ( rho_ice * r_ice**3. ) dqrcldfreez = MAX(dqrcldfreez, - raincld(i) / dhum_to_dflux(i)) dqrcldfreez = MAX(dqrcldfreez, dqrfreez_max) dqrtotfreez_step1 = dqrcldfreez !--Add tendencies !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision) qice(i) = qice(i) + dqifreez raincld(i) = MAX(0., raincld(i) + dqrcldfreez * dhum_to_dflux(i)) snowcld(i) = MAX(0., snowcld(i) - dqrcldfreez * dhum_to_dflux(i) - dqifreez * dhum_to_dflux(i)) temp(i) = temp(i) - dqrtotfreez_step1 * RLMLT / RCPD & / ( 1. + RVTMP2 * qtot(i) ) ENDIF !-- Second step immersion freezing of rain drops !-- with a homemade timeconstant depending on temperature dqrclrfreez = 0. dqrcldfreez = 0. dqrtotfreez_step2 = 0. !--Computed according to !--Cpdry * Delta T * (1 + (Cpvap/Cpdry - 1) * qtot) = Lfusion * Delta q dqrfreez_max = MIN(0., ( temp(i) - RTT ) / RLMLT * RCPD & * ( 1. + RVTMP2 * qtot(i) )) tau_freez = 1. / ( beta_freez & * EXP( - alpha_freez * ( temp(i) - temp_nowater ) / ( RTT - temp_nowater ) ) ) !--In clear air IF ( rainclr(i) .GT. 0. ) THEN !--Exact solution of dqrain/dt = -qrain/tau_freez dqrclrfreez = rainclr(i) / dhum_to_dflux(i) * ( EXP( - dtime / tau_freez ) - 1. ) ENDIF !--In cloudy air IF ( raincld(i) .GT. 0. ) THEN !--Exact solution of dqrain/dt = -qrain/tau_freez dqrcldfreez = raincld(i) / dhum_to_dflux(i) * ( EXP( - dtime / tau_freez ) - 1. ) ENDIF !--temperature barrier step 2 !--It is activated if the total is LOWER than the max !--because everything is negative dqrtotfreez_step2 = dqrclrfreez + dqrcldfreez IF ( dqrtotfreez_step2 .LT. dqrfreez_max ) THEN !--We redistribute the max freezed rain keeping !--the clear/cloud partition of the freezing rain dqrclrfreez = dqrfreez_max * dqrclrfreez / dqrtotfreez_step2 dqrcldfreez = dqrfreez_max * dqrcldfreez / dqrtotfreez_step2 dqrtotfreez_step2 = dqrfreez_max ENDIF !--Add tendencies !--The MAX is needed because in some cases, the flux can be slightly negative (numerical precision) rainclr(i) = MAX(0., rainclr(i) + dqrclrfreez * dhum_to_dflux(i)) raincld(i) = MAX(0., raincld(i) + dqrcldfreez * dhum_to_dflux(i)) snowclr(i) = MAX(0., snowclr(i) - dqrclrfreez * dhum_to_dflux(i)) snowcld(i) = MAX(0., snowcld(i) - dqrcldfreez * dhum_to_dflux(i)) !--Temperature adjustment with the uptake of latent !--heat because of freezing temp(i) = temp(i) - dqrtotfreez_step2 * RLMLT / RCPD & / ( 1. + RVTMP2 * qtot(i) ) !--Diagnostic tendencies dqrtotfreez = dqrtotfreez_step1 + dqrtotfreez_step2 dqrfreez(i) = dqrtotfreez / dtime dqsfreez(i) = -(dqrtotfreez + dqifreez) / dtime ENDIF !--If the local flux of rain+snow in clear air is lower than min_precip, !--we reduce the precipiration fraction in the clear air so that the new !--local flux of rain+snow is equal to min_precip. !--we apply the minimum only on the clear-sky fraction because the cloudy precip fraction !--already decreases out of clouds !--Here, rain+snow is the gridbox-mean flux of precip. !--Therefore, (rain+snow)/precipfrac is the local flux of precip. !--If the local flux of precip is lower than min_precip, i.e., !-- (rain+snow)/precipfrac < min_precip , i.e., !-- (rain+snow)/min_precip < precipfrac , then we want to reduce !--the precip fraction to the equality, i.e., precipfrac = (rain+snow)/min_precip. !--Note that this is physically different than what is proposed in LTP thesis. !--min_precip is either equal to rain_int_min or calculated as a very small fraction !--of the minimum precip flux estimated as the flux associated with the !--autoconversion threshold mass content !min_precip=1.e-6*(pplay(i)/RD/temp(i))*MIN(rain_fallspeed_clr*cld_lc_lsc,snow_fallspeed_clr*cld_lc_lsc_snow) min_precip=rain_int_min precipfracclr(i) = MIN( precipfracclr(i), ( rainclr(i) + snowclr(i) ) / min_precip ) !--Calculate outputs rain(i) = rainclr(i) + raincld(i) snow(i) = snowclr(i) + snowcld(i) !--Diagnostics !--BEWARE this is indeed a diagnostic: this is an estimation from !--the value of the flux at the bottom interface of the mesh and !--and assuming an upstream numerical calculation !--If ok_radocond_snow is TRUE, then the diagnostic qsnowdiag is !--used for computing the total ice water content in the mesh !--for radiation only qraindiag(i) = ( rainclr(i) / ( pplay(i) / RD / temp(i) ) / rain_fallspeed_clr & + raincld(i) / ( pplay(i) / RD / temp(i) ) / rain_fallspeed_cld ) qsnowdiag(i) = ( snowclr(i) / ( pplay(i) / RD / temp(i) ) / snow_fallspeed_clr & + snowcld(i) / ( pplay(i) / RD / temp(i) ) / snow_fallspeed_cld ) ENDDO ! loop on klon END SUBROUTINE poprecip_postcld END MODULE lmdz_lscp_precip