MODULE LSCP_mod IMPLICIT NONE CONTAINS !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE LSCP(dtime,missing_val, & paprs,pplay,t,q,ptconv,ratqs, & d_t, d_q, d_ql, d_qi, rneb, rneb_seri, & radliq, radicefrac, rain, snow, & pfrac_impa, pfrac_nucl, pfrac_1nucl, & frac_impa, frac_nucl, beta, & prfl, psfl, rhcl, zqta, fraca, & ztv, zpspsk, ztla, zthl, iflag_cld_th, & iflag_ice_thermo, ok_ice_sursat) !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! Authors: Z.X. Li (LMD), J-L Dufresne (LMD), C. Rio (LMD), J-Y Grandpeix (LMD) ! A. JAM (LMD), J-B Madeleine (LMD), E. Vignon (LMD), L. Touzze-Peiffert (LMD) !-------------------------------------------------------------------------------- ! Date: 01/2021 !-------------------------------------------------------------------------------- ! Aim: Large Scale Clouds and Precipitation (LSCP) ! ! This code is a new version of the fisrtilp.F90 routine, which is itself a ! merge of 'first' (superrsaturation physics, P. LeVan K. Laval) ! and 'ilp' (il pleut, L. Li) ! ! Compared to fisrtilp, LSCP ! -> assumes thermcep = .TRUE. all the time (fisrtilp inconsistent when .FALSE.) ! -> consider always precipitation thermalisation (fl_cor_ebil>0) ! -> option iflag_fisrtilp_qsat<0 no longer possible (qsat does not evolve with T) ! -> option "oldbug" by JYG has been removed ! -> iflag_t_glace >0 always ! -> the 'all or nothing' cloud approach is no longer available (cpartiel=T always) ! -> rectangular distribution from L. Li no longer available ! -> We always account for the Wegener-Findeisen-Bergeron process (iflag_bergeron = 2 in fisrt) !-------------------------------------------------------------------------------- ! References: ! ! - Hourdin et al. 2013, Clim Dyn, doi:10.1007/s00382-012-1343-y ! - Jam et al. 2013, Boundary-Layer Meteorol, doi:10.1007/s10546-012-9789-3 ! - Madeleine et al. 2020, JAMES, doi:10.1029/2020MS002046 ! ------------------------------------------------------------------------------- ! Code structure: ! ! P0> Thermalization of the precipitation coming from the overlying layer ! P1> Evaporation of the precipitation (falling from the k+1 level) ! P2> Cloud formation (at the k level) ! P2.A.0> Cloud properties calculation from a rectangular pdf ! P2.A.1> With the new PDFs, calculation of cloud properties using the inital ! values of T and Q ! P2.A.2> Coupling between condensed water and temperature ! P2.A.3> Calculation of final quantities associated with cloud formation ! P2.B> 'All or nothing' cloud ! P2.C> Release of Latent heat after cloud formation ! P3> Autoconversion to precipitation (k-level) ! P4> Wet scavenging !------------------------------------------------------------------------------ ! Some preliminary comments (JBM) : ! ! The cloud water that the radiation scheme sees is not the same that the cloud ! water used in the physics and the dynamics ! ! During the autoconversion to precipitation (P3 step), radliq (cloud water used ! by the radiation scheme) is calculated as an average of the water that remains ! in the cloud during the precipitation and not the water remaining at the end ! of the time step. The latter is used in the rest of the physics and advected ! by the dynamics. ! ! In summary: ! ! Radiation: ! xflwc(newmicro)+xfiwc(newmicro) = ! cldliq(physiq)=radliq(fisrt)=lwcon(nc)+iwcon(nc) ! ! Physics/Dynamics: ! ql_seri(physiq)+qs_seri(physiq)=ocond(nc) ! ! Notetheless, be aware of: ! ! radliq(fisrt) .NE. ocond(nc) ! i.e.: ! lwcon(nc)+iwcon(nc) .NE. ocond(nc) ! (which is not trivial) !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ! USE dimphy USE print_control_mod, ONLY: prt_level, lunout USE cloudth_mod USE ioipsl_getin_p_mod, ONLY : getin_p USE phys_local_var_mod, ONLY: ql_seri,qs_seri USE phys_local_var_mod, ONLY: rneblsvol USE lscp_tools_mod, ONLY : CALC_QSAT_ECMWF, ICEFRAC_LSCP, CALC_GAMMASAT, FALLICE_VELOCITY USE ice_sursat_mod !--ice supersaturation USE phys_local_var_mod, ONLY: zqsats, zqsatl USE phys_local_var_mod, ONLY: qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss USE phys_local_var_mod, ONLY: Tcontr, qcontr, qcontr2, fcontrN, fcontrP IMPLICIT NONE !=============================================================================== ! VARIABLE DECLARATION !=============================================================================== include "YOMCST.h" include "YOETHF.h" include "FCTTRE.h" include "fisrtilp.h" include "nuage.h" ! INPUT VARIABLES: !----------------- REAL, INTENT(IN) :: dtime ! time step [s] REAL, INTENT(IN) :: missing_val ! missing value for output REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! inter-layer pressure [Pa] REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! mid-layer pressure [Pa] REAL, DIMENSION(klon,klev), INTENT(IN) :: t ! temperature (K) REAL, DIMENSION(klon,klev), INTENT(IN) :: q ! specific humidity [kg/kg] INTEGER, INTENT(IN) :: iflag_cld_th ! flag that determines the distribution of convective clouds INTEGER, INTENT(IN) :: iflag_ice_thermo! flag to activate the ice thermodynamics ! CR: if iflag_ice_thermo=2, only convection is active LOGICAL, INTENT(IN) :: ok_ice_sursat ! flag to determine if ice sursaturation is activated LOGICAL, DIMENSION(klon,klev), INTENT(IN) :: ptconv ! grid points where deep convection scheme is active !Inputs associated with thermal plumes REAL, DIMENSION(klon,klev), INTENT(IN) :: ztv ! virtual potential temperature [K] REAL, DIMENSION(klon,klev), INTENT(IN) :: zqta ! specific humidity within thermals [kg/kg] REAL, DIMENSION(klon,klev), INTENT(IN) :: fraca ! fraction of thermals within the mesh [-] REAL, DIMENSION(klon,klev), INTENT(IN) :: zpspsk ! exner potential (p/100000)**(R/cp) REAL, DIMENSION(klon,klev), INTENT(IN) :: ztla ! liquid temperature within thermals [K] REAL, DIMENSION(klon,klev), INTENT(INOUT) :: zthl ! liquid potential temperature [K] ! INPUT/OUTPUT variables !------------------------ REAL, DIMENSION(klon,klev), INTENT(INOUT):: ratqs ! function of pressure that sets the large-scale ! Input sursaturation en glace REAL, DIMENSION(klon,klev), INTENT(INOUT):: rneb_seri ! fraction nuageuse en memoire ! OUTPUT variables !----------------- REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_t ! temperature increment [K] REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_q ! specific humidity increment [kg/kg] REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_ql ! liquid water increment [kg/kg] REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_qi ! cloud ice mass increment [kg/kg] REAL, DIMENSION(klon,klev), INTENT(OUT) :: rneb ! cloud fraction [-] REAL, DIMENSION(klon,klev), INTENT(OUT) :: radliq ! condensed water used in the radiation scheme [kg/kg] REAL, DIMENSION(klon,klev), INTENT(OUT) :: radicefrac ! ice fraction of condensed water for radiation scheme REAL, DIMENSION(klon,klev), INTENT(OUT) :: rhcl ! clear-sky relative humidity [-] REAL, DIMENSION(klon), INTENT(OUT) :: rain ! large-scale rainfall [kg/s/m2] REAL, DIMENSION(klon), INTENT(OUT) :: snow ! large-scale snowfall [kg/s/m2] REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: prfl ! large-scale rainfall flux in the column [kg/s/m2] REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: psfl ! large-scale snowfall flux in the column [kg/s/m2] ! cofficients of scavenging fraction (for off-line) REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfrac_nucl ! scavenging fraction due tu nucleation [-] REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfrac_1nucl ! scavenging fraction due tu nucleation with a -1 factor [-] REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfrac_impa ! scavening fraction due to impaction [-] ! fraction of aerosol scavenging through impaction and nucleation (for on-line) REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_impa ! scavenging fraction due tu impaction [-] REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_nucl ! scavenging fraction due tu nucleation [-] ! PROGRAM PARAMETERS: !-------------------- REAL, SAVE :: seuil_neb=0.001 ! cloud fraction threshold: a cloud really exists when exceeded !$OMP THREADPRIVATE(seuil_neb) INTEGER, SAVE :: ninter=5 ! number of iterations to calculate autoconversion to precipitation !$OMP THREADPRIVATE(ninter) INTEGER,SAVE :: iflag_evap_prec=1 ! precipitation evaporation flag. 0: nothing, 1: "old way", ! 2: Max cloud fraction above to calculate the max of reevaporation ! 4: LTP'method i.e. evaporation in the clear-sky fraction of the mesh only !$OMP THREADPRIVATE(iflag_evap_prec) REAL t_coup ! temperature threshold which determines the phase PARAMETER (t_coup=234.0) ! for which the saturation vapor pressure is calculated REAL DDT0 ! iteration parameter PARAMETER (DDT0=.01) REAL ztfondue ! parameter to calculate melting fraction of precipitation PARAMETER (ztfondue=278.15) REAL, SAVE :: rain_int_min=0.001 ! Minimum local rain intensity [mm/s] before the decrease in associated precipitation fraction !$OMP THREADPRIVATE(rain_int_min) ! LOCAL VARIABLES: !---------------- REAL qsl(klon), qsi(klon) REAL zct, zcl REAL ctot(klon,klev) REAL ctot_vol(klon,klev) REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5 REAL zdqsdT_raw(klon) REAL gammasat(klon),dgammasatdt(klon) ! coefficient to make cold condensation at the correct RH and derivative wrt T REAL Tbef(klon),qlbef(klon),DT(klon),num,denom REAL cste REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon) REAL Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon) REAL erf REAL qcloud(klon), icefrac_mpc(klon,klev), qincloud_mpc(klon) REAL zrfl(klon), zrfln(klon), zqev, zqevt REAL zifl(klon), zifln(klon), ziflprev(klon),zqev0,zqevi, zqevti REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon) REAL zoliql(klon), zoliqi(klon) REAL zt(klon) REAL zdz(klon),zrho(klon,klev),iwc(klon) REAL zchau,zfroi,zfice(klon),zneb(klon),znebprecip(klon) REAL zmelt,zrain,zsnow,zprecip REAL dzfice(klon) REAL zsolid REAL qtot(klon), qzero(klon) REAL dqsl(klon),dqsi(klon) REAL smallestreal ! Variables for Bergeron process REAL zcp, coef1, DeltaT, Deltaq, Deltaqprecl REAL zqpreci(klon), zqprecl(klon) ! Variables precipitation energy conservation REAL zmqc(klon) ! Variables for tracers ! temporary: alpha parameter for scavenging ! 4 possible scavenging processes REAL,SAVE :: a_tr_sca(4) !$OMP THREADPRIVATE(a_tr_sca) REAL zalpha_tr REAL zfrac_lessi REAL zprec_cond(klon) REAL beta(klon,klev) ! conversion rate of condensed water REAL zmair(klon), zcpair, zcpeau REAL zlh_solid(klon), zm_solid ! for liquid -> solid conversion REAL zrflclr(klon), zrflcld(klon) REAL d_zrfl_clr_cld(klon), d_zifl_clr_cld(klon) REAL d_zrfl_cld_clr(klon), d_zifl_cld_clr(klon) REAL ziflclr(klon), ziflcld(klon) REAL znebprecipclr(klon), znebprecipcld(klon) REAL tot_zneb(klon), tot_znebn(klon), d_tot_zneb(klon) REAL d_znebprecip_clr_cld(klon), d_znebprecip_cld_clr(klon) REAL velo(klon,klev), vr REAL qlmpc, qimpc, rnebmpc REAL radliqi(klon,klev), radliql(klon,klev) INTEGER i, k, n, kk INTEGER n_i(klon), iter INTEGER ncoreczq INTEGER mpc_bl_points(klon,klev) INTEGER,SAVE :: itap=0 !$OMP THREADPRIVATE(itap) LOGICAL lognormale(klon) LOGICAL convergence(klon) LOGICAL,SAVE :: appel1er !$OMP THREADPRIVATE(appel1er) CHARACTER (len = 20) :: modname = 'lscp' CHARACTER (len = 80) :: abort_message DATA appel1er /.TRUE./ !=============================================================================== ! INITIALISATION !=============================================================================== ! Few initial checks IF (iflag_t_glace.EQ.0) THEN abort_message = 'lscp cannot be used if iflag_t_glace=0' CALL abort_physic(modname,abort_message,1) ENDIF IF (.NOT.((iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3))) THEN abort_message = 'lscp cannot be used without ice thermodynamics' CALL abort_physic(modname,abort_message,1) ENDIF IF (.NOT.thermcep) THEN abort_message = 'lscp cannot be used when thermcep=false' CALL abort_physic(modname,abort_message,1) ENDIF IF (iflag_fisrtilp_qsat .LT. 0) THEN abort_message = 'lscp cannot be used with iflag_fisrtilp<0' CALL abort_physic(modname,abort_message,1) ENDIF ! Few initialisations itap=itap+1 znebprecip(:)=0.0 ctot_vol(1:klon,1:klev)=0.0 rneblsvol(1:klon,1:klev)=0.0 smallestreal=1.e-9 znebprecipclr(:)=0.0 znebprecipcld(:)=0.0 mpc_bl_points(:,:)=0 IF (prt_level>9) WRITE(lunout,*) 'NUAGES4 A. JAM' IF (appel1er) THEN CALL getin_p('ninter',ninter) CALL getin_p('iflag_evap_prec',iflag_evap_prec) CALL getin_p('seuil_neb',seuil_neb) CALL getin_p('rain_int_min',rain_int_min) WRITE(lunout,*) 'lscp, ninter:', ninter WRITE(lunout,*) 'lscp, iflag_evap_prec:', iflag_evap_prec WRITE(lunout,*) 'lscp, seuil_neb:', seuil_neb WRITE(lunout,*) 'lscp, rain_int_min:', rain_int_min ! check for precipitation sub-time steps IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN WRITE(lunout,*) 'lscp: it is not expected, see Z.X.Li', dtime WRITE(lunout,*) 'I would prefer a 6 min sub-timestep' ENDIF !AA Temporary initialisation a_tr_sca(1) = -0.5 a_tr_sca(2) = -0.5 a_tr_sca(3) = -0.5 a_tr_sca(4) = -0.5 !AA Set scavenged fractions to 1 DO k = 1, klev DO i = 1, klon pfrac_nucl(i,k)=1.0 pfrac_1nucl(i,k)=1.0 pfrac_impa(i,k)=1.0 beta(i,k)=0.0 ENDDO ENDDO IF (ok_ice_sursat) CALL ice_sursat_init() appel1er = .FALSE. ENDIF !(appel1er) ! AA for 'safety' reasons zalpha_tr = 0. zfrac_lessi = 0. ! Initialisation of output variables (JYG): prfl(:,:) = 0.0 psfl(:,:) = 0.0 d_t(:,:) = 0.0 d_q(:,:) = 0.0 d_ql(:,:) = 0.0 d_qi(:,:) = 0.0 rneb(:,:) = 0.0 radliq(:,:) = 0.0 radicefrac(:,:) = 0.0 frac_nucl(:,:) = 1.0 frac_impa(:,:) = 1.0 rain(:) = 0.0 snow(:) = 0.0 zoliq(:)=0.0 zfice(:)=0.0 dzfice(:)=0.0 zqprecl(:)=0.0 zqpreci(:)=0.0 zrfl(:) = 0.0 zifl(:) = 0.0 ziflprev(:)=0.0 zneb(:) = seuil_neb zrflclr(:) = 0.0 ziflclr(:) = 0.0 zrflcld(:) = 0.0 ziflcld(:) = 0.0 tot_zneb(:) = 0.0 tot_znebn(:) = 0.0 d_tot_zneb(:) = 0.0 qzero(:) = 0.0 !--ice sursaturation gamma_ss(:,:) = 1.0 qss(:,:) = 0.0 rnebss(:,:) = 0.0 Tcontr(:,:) = missing_val qcontr(:,:) = missing_val qcontr2(:,:) = missing_val fcontrN(:,:) = 0.0 fcontrP(:,:) = 0.0 !=============================================================================== ! BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM !=============================================================================== ncoreczq=0 DO k = klev, 1, -1 ! Initialisation temperature and specific humidity DO i = 1, klon zt(i)=t(i,k) zq(i)=q(i,k) ENDDO ! -------------------------------------------------------------------- ! P0> Thermalization of precipitation falling from the overlying layer ! -------------------------------------------------------------------- ! Computes air temperature variation due to latent heat 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 (k.LE.klevm1) THEN DO i = 1, klon zmair(i)=(paprs(i,k)-paprs(i,k+1))/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(i) ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) & / (zcpair + zmqc(i)*zcpeau) ENDDO ELSE DO i = 1, klon zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG zmqc(i) = 0. ENDDO ENDIF ! -------------------------------------------------------------------- ! P1> Precipitation evaporation/sublimation ! -------------------------------------------------------------------- ! A part of the precipitation coming from above is evaporated/sublimated. ! -------------------------------------------------------------------- IF (iflag_evap_prec.GE.1) THEN ! Calculation of saturation specific humidity ! depending on temperature: CALL CALC_QSAT_ECMWF(klon,zt(:),qzero(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:)) ! wrt liquid water CALL CALC_QSAT_ECMWF(klon,zt(:),qzero(:),pplay(:,k),RTT,1,.false.,qsl(:),dqsl(:)) ! wrt ice CALL CALC_QSAT_ECMWF(klon,zt(:),qzero(:),pplay(:,k),RTT,2,.false.,qsi(:),dqsi(:)) DO i = 1, klon ! if precipitation IF (zrfl(i)+zifl(i).GT.0.) THEN ! LTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec=4). IF (iflag_evap_prec.EQ.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.EQ.4) THEN ! Max evaporation not to saturate the whole mesh 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))) & *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG ELSE IF (iflag_evap_prec.EQ.4) THEN zqevt = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsl(i)) & *SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) & *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG ELSE zqevt = 1.*coef_eva*(1.0-zq(i)/qsl(i))*SQRT(zrfl(i)) & *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG ENDIF zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) & *RG*dtime/(paprs(i,k)-paprs(i,k+1)) ! sublimation of the solid precipitation coming from above IF (iflag_evap_prec.EQ.3) THEN zqevti = znebprecip(i)*coef_eva_i*(1.0-zq(i)/qsi(i)) & *SQRT(zifl(i)/max(1.e-4,znebprecip(i))) & *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG ELSE IF (iflag_evap_prec.EQ.4) THEN zqevti = znebprecipclr(i)*coef_eva_i*(1.0-zq(i)/qsi(i)) & *SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) & *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG ELSE zqevti = 1.*coef_eva_i*(1.0-zq(i)/qsi(i))*SQRT(zifl(i)) & *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG ENDIF zqevti = MAX(0.0,MIN(zqevti,zifl(i))) & *RG*dtime/(paprs(i,k)-paprs(i,k+1)) ! 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 ! New solid and liquid precipitation fluxes zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) & /RG/dtime) zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) & /RG/dtime) ! vapor, temperature, precip fluxes update zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) & * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime ! precip thermalization zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) & * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime zt(i) = zt(i) + (zrfln(i)-zrfl(i)) & * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime & * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) & + (zifln(i)-zifl(i)) & * (RG/(paprs(i,k)-paprs(i,k+1)))*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) IF (iflag_evap_prec.EQ.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 ! CR Be careful: ice melting has been moved zmelt = ((zt(i)-273.15)/(ztfondue-273.15)) ! JYG ! precip fraction that is melted zmelt = MIN(MAX(zmelt,0.),1.) ! Ice melting IF (iflag_evap_prec.EQ.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 ! Latent heat of melting with precipitation thermalization zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) & *RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) IF (iflag_evap_prec.EQ.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 ENDIF ! (iflag_evap_prec>=1) ! -------------------------------------------------------------------- ! End precip evaporation ! -------------------------------------------------------------------- ! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter !------------------------------------------------------- qtot(:)=zq(:)+zmqc(:) CALL CALC_QSAT_ECMWF(klon,zt(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:)) DO i = 1, klon zdelta = MAX(0.,SIGN(1.,RTT-zt(i))) zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta) IF (zq(i) .LT. 1.e-15) THEN ncoreczq=ncoreczq+1 zq(i)=1.e-15 ENDIF ENDDO ! -------------------------------------------------------------------- ! P2> Cloud formation !--------------------------------------------------------------------- ! ! Unlike fisrtilp, we always assume a 'fractional cloud' approach ! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution ! is prescribed and depends on large scale variables and boundary layer ! properties) ! The decrease in condensed part due tu latent heating is taken into ! account ! ------------------------------------------------------------------- ! P2.1> With the PDFs (log-normal, bigaussian) ! cloud propertues calculation with the initial values of t and q ! ---------------------------------------------------------------- ! initialise gammasat and qincloud_mpc gammasat(:)=1. qincloud_mpc(:)=0. IF (iflag_cld_th.GE.5) THEN IF (iflag_mpc_bl .LT. 1) THEN IF (iflag_cloudth_vert.LE.2) THEN CALL cloudth(klon,klev,k,ztv, & zq,zqta,fraca, & qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, & ratqs,zqs,t) ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN CALL cloudth_v3(klon,klev,k,ztv, & zq,zqta,fraca, & qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, & ratqs,zqs,t) !Jean Jouhaud's version, Decembre 2018 !------------------------------------- ELSEIF (iflag_cloudth_vert.EQ.6) THEN CALL cloudth_v6(klon,klev,k,ztv, & zq,zqta,fraca, & qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, & ratqs,zqs,t) ENDIF ELSE ! cloudth_v3 + cold microphysical considerations ! + stationary mixed-phase cloud model CALL cloudth_mpc(klon,klev,k,mpc_bl_points, & t,ztv,zq,zqta,fraca, zpspsk, & paprs,pplay,ztla,zthl,ratqs,zqs,psfl, & qcloud,qincloud_mpc,icefrac_mpc,ctot,ctot_vol) ENDIF ! iflag_mpc_bl DO i=1,klon rneb(i,k)=ctot(i,k) rneblsvol(i,k)=ctot_vol(i,k) zqn(i)=qcloud(i) ENDDO ENDIF IF (iflag_cld_th .LE. 4) THEN ! lognormal lognormale = .TRUE. ELSEIF (iflag_cld_th .GE. 6) THEN ! lognormal distribution when no thermals lognormale = fraca(:,k) < 1e-10 ELSE ! When iflag_cld_th=5, we always assume ! bi-gaussian distribution lognormale = .FALSE. ENDIF DT(:) = 0. n_i(:)=0 Tbef(:)=zt(:) qlbef(:)=0. ! Treatment of non-boundary layer clouds (lognormale) ! condensation with qsat(T) variation (adaptation) ! Iterative Loop to converge towards qsat DO iter=1,iflag_fisrtilp_qsat+1 DO i=1,klon ! convergence = .true. until when convergence is not satisfied convergence(i)=ABS(DT(i)).GT.DDT0 IF ((convergence(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN ! if not convergence: ! P2.2.1> cloud fraction and condensed water mass calculation ! Calculated variables: ! rneb : cloud fraction ! zqn : total water within the cloud ! zcond: mean condensed water within the mesh ! rhcl: clear-sky relative humidity !--------------------------------------------------------------- ! new temperature: Tbef(i)=Tbef(i)+DT(i) ! Rneb, qzn and zcond for lognormal PDFs qtot(i)=zq(i)+zmqc(i) ENDIF ENDDO ! Calculation of saturation specific humidity and ce fraction CALL calc_qsat_ecmwf(klon,Tbef(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:)) CALL calc_gammasat(klon,Tbef(:),qtot(:),pplay(:,k),gammasat(:),dgammasatdt(:)) ! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:) CALL icefrac_lscp(klon, zt(:),pplay(:,k)/paprs(:,1),zfice(:),dzfice(:)) DO i=1,klon IF ((convergence(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN zpdf_sig(i)=ratqs(i,k)*zq(i) zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2)) zpdf_delta(i)=log(zq(i)/(gammasat(i)*zqs(i))) zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.)) zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.)) zpdf_e1(i)=zpdf_a(i)-zpdf_b(i) zpdf_e1(i)=sign(min(ABS(zpdf_e1(i)),5.),zpdf_e1(i)) zpdf_e1(i)=1.-erf(zpdf_e1(i)) zpdf_e2(i)=zpdf_a(i)+zpdf_b(i) zpdf_e2(i)=sign(min(ABS(zpdf_e2(i)),5.),zpdf_e2(i)) zpdf_e2(i)=1.-erf(zpdf_e2(i)) !--ice sursaturation by Audran IF ((.NOT.ok_ice_sursat).OR.(Tbef(i).GT.t_glace_min)) THEN IF (zpdf_e1(i).LT.1.e-10) THEN rneb(i,k)=0. zqn(i)=gammasat(i)*zqs(i) ELSE rneb(i,k)=0.5*zpdf_e1(i) zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i) ENDIF rnebss(i,k)=0.0 !--added by OB (needed because of the convergence loop on time) fcontrN(i,k)=0.0 !--idem fcontrP(i,k)=0.0 !--idem qss(i,k)=0.0 !--idem ELSE !------------------------------------ ! ICE SUPERSATURATION !------------------------------------ CALL ice_sursat(pplay(i,k), paprs(i,k)-paprs(i,k+1), dtime, i, k, t(i,k), zq(i), & gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k), & rneb(i,k), zqn(i), rnebss(i,k), qss(i,k), & Tcontr(i,k), qcontr(i,k), qcontr2(i,k), fcontrN(i,k), fcontrP(i,k) ) ENDIF ! ((flag_ice_sursat.eq.0).or.(Tbef(i).gt.t_glace_min)) ! If vertical heterogeneity, change fraction by volume as well IF (iflag_cloudth_vert.GE.3) THEN ctot_vol(i,k)=rneb(i,k) rneblsvol(i,k)=ctot_vol(i,k) ENDIF ! P2.2.2> Approximative calculation of temperature variation DT ! due to condensation. ! Calculated variables: ! dT : temperature change due to condensation !--------------------------------------------------------------- IF (zfice(i).LT.1) THEN cste=RLVTT ELSE cste=RLSTT ENDIF qlbef(i)=max(0.,zqn(i)-zqs(i)) num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT & +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i) denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) & -(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k) & *qlbef(i)*dzfice(i) DT(i)=num/denom n_i(i)=n_i(i)+1 ENDIF ! end convergence ENDDO ! end loop on i ENDDO ! iter=1,iflag_fisrtilp_qsat+1 ! P2.3> Final quantities calculation ! Calculated variables: ! rneb : cloud fraction ! zcond: mean condensed water in the mesh ! zqn : mean water vapor in the mesh ! zt : temperature ! rhcl : clear-sky relative humidity ! ---------------------------------------------------------------- ! Water vapor update, Phase determination and subsequent latent heat exchange ! Partition function in stratiform clouds (will be overwritten in boundary-layer MPCs) CALL icefrac_lscp(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:), dzfice(:)) DO i=1, klon IF (mpc_bl_points(i,k) .GT. 0) THEN zcond(i) = MAX(0.0,qincloud_mpc(i))*rneb(i,k) ! following line is very strange and probably wrong rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i) ! water vapor update and partition function if thermals zq(i) = zq(i) - zcond(i) zfice(i)=icefrac_mpc(i,k) ELSE ! Checks on rneb, rhcl and zqn IF (rneb(i,k) .LE. 0.0) THEN zqn(i) = 0.0 rneb(i,k) = 0.0 zcond(i) = 0.0 rhcl(i,k)=zq(i)/zqs(i) ELSE IF (rneb(i,k) .GE. 1.0) THEN zqn(i) = zq(i) rneb(i,k) = 1.0 zcond(i) = MAX(0.0,zqn(i)-gammasat(i)*zqs(i)) rhcl(i,k)=1.0 ELSE zcond(i) = MAX(0.0,zqn(i)-gammasat(i)*zqs(i))*rneb(i,k) ! following line is very strange and probably wrong: rhcl(i,k)=(zqs(i)+zq(i))/2./zqs(i) ENDIF ! water vapor update zq(i) = zq(i) - zcond(i) ENDIF ! temperature update due to phase change zt(i) = zt(i) + (1.-zfice(i))*zcond(i) & & * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) & +zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) ENDDO ! If vertical heterogeneity, change volume fraction IF (iflag_cloudth_vert .GE. 3) THEN ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.) rneblsvol(1:klon,k)=ctot_vol(1:klon,k) ENDIF ! ---------------------------------------------------------------- ! P3> Precipitation formation ! ---------------------------------------------------------------- ! LTP: IF (iflag_evap_prec==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,k),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,k),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,k) - 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) ENDDO ENDIF ! Initialisation of zoliq and radliq variables DO i = 1, klon zoliq(i) = zcond(i) zoliqi(i)= zoliq(i)*zfice(i) zoliql(i)= zoliq(i)*(1.-zfice(i)) zrho(i,k) = pplay(i,k) / zt(i) / RD zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG) iwc(i) = 0. zneb(i) = MAX(rneb(i,k),seuil_neb) radliq(i,k) = zoliq(i)/REAL(ninter+1) radicefrac(i,k) = zfice(i) radliqi(i,k)=zoliq(i)*zfice(i)/REAL(ninter+1) radliql(i,k)=zoliq(i)*(1.-zfice(i))/REAL(ninter+1) ENDDO ! ------------------------------------------------------------------------------- ! Autoconversion ! ------------------------------------------------------------------------------- DO n = 1, ninter DO i=1,klon IF (rneb(i,k).GT.0.0) THEN iwc(i) = zrho(i,k) * zoliqi(i) / zneb(i) ! in-cloud ice condensate content ENDIF ENDDO CALL fallice_velocity(klon,iwc(:),zt(:),zrho(:,k),pplay(:,k),ptconv(:,k),velo(:,k)) DO i = 1, klon IF (rneb(i,k).GT.0.0) THEN ! Initialization of zrain, zsnow and zprecip: zrain=0. zsnow=0. zprecip=0. IF (zneb(i).EQ.seuil_neb) THEN zprecip = 0.0 zsnow = 0.0 zrain= 0.0 ELSE ! water quantity to remove: zchau (Sundqvist, 1978) ! same thing for the ice: zfroi (Zender & Kiehl, 1997) IF (ptconv(i,k)) THEN ! if convective point zcl=cld_lc_con zct=1./cld_tau_con zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i)*velo(i,k)*zfice(i) ELSE zcl=cld_lc_lsc zct=1./cld_tau_lsc zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) & ! dqice/dt=1/rho*d(rho*wice*qice)/dz *velo(i,k) * zfice(i) 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 zchau = zct *dtime/REAL(ninter) * zoliq(i) & *(1.0-EXP(-(zoliq(i)/ctot_vol(i,k)/zcl)**2)) *(1.-zfice(i)) ELSE zchau = zct *dtime/REAL(ninter) * zoliq(i) & *(1.0-EXP(-(zoliq(i)/zneb(i)/zcl)**2)) *(1.-zfice(i)) ! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2) ENDIF zrain = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i))) zsnow = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i)) zprecip = MAX(zrain + zsnow,0.0) ENDIF 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) radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1) radliql(i,k) = radliql(i,k) + zoliql(i)/REAL(ninter+1) radliqi(i,k) = radliqi(i,k) + zoliqi(i)/REAL(ninter+1) ENDIF ! rneb >0 ENDDO ! i = 1,klon ENDDO ! n = 1,niter ! Precipitation flux calculation DO i = 1, klon ziflprev(i)=zifl(i) IF (rneb(i,k) .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 zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i) zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i)) zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) coef1 = rneb(i,k)*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,k)*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.) IF (iflag_evap_prec.EQ.4) THEN zrflcld(i) = zrflcld(i)+zqprecl(i) & *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) ziflcld(i) = ziflcld(i)+ zqpreci(i) & *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) znebprecipcld(i) = rneb(i,k) zrfl(i) = zrflcld(i) + zrflclr(i) zifl(i) = ziflcld(i) + ziflclr(i) ELSE zrfl(i) = zrfl(i)+ zqprecl(i) & *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) zifl(i) = zifl(i)+ zqpreci(i) & *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) ENDIF 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_pre=4 IF (iflag_evap_prec.EQ.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 ! End of precipitation formation ! -------------------------------- ! Outputs: ! Precipitation fluxes at layer interfaces ! and temperature and water species tendencies DO i = 1, klon psfl(i,k)=zifl(i) prfl(i,k)=zrfl(i) d_ql(i,k) = (1-zfice(i))*zoliq(i) d_qi(i,k) = zfice(i)*zoliq(i) d_q(i,k) = zq(i) - q(i,k) d_t(i,k) = zt(i) - t(i,k) ENDDO ! Calculation of the concentration of condensates seen by the radiation scheme ! for liquid phase, we take radliql ! for ice phase, we take radliqi if we neglect snowfall, otherwise (ok_radliq_snow=true) ! we recaulate radliqi to account for contributions from the precipitation flux IF ((ok_radliq_snow) .AND. (k .LT. klev)) THEN ! for the solid phase (crystals + snowflakes) ! we recalculate radliqi 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 radliqi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)+ziflprev(i)/zrho(i,k+1)/velo(i,k+1) ELSE radliqi(i,k)=zoliq(i)*zfice(i)+zqpreci(i) ENDIF radliq(i,k)=radliql(i,k)+radliqi(i,k) ENDDO ENDIF ! caculate the percentage of ice in "radliq" so cloud+precip seen by the radiation scheme DO i=1,klon IF (radliq(i,k) .GT. 0.) THEN radicefrac(i,k)=MIN(MAX(radliqi(i,k)/radliq(i,k),0.),1.) ENDIF ENDDO ! ---------------------------------------------------------------- ! P4> Wet scavenging ! ---------------------------------------------------------------- !Scavenging through nucleation in the layer DO i = 1,klon IF(zcond(i).GT.zoliq(i)+1.e-10) THEN beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime ELSE beta(i,k) = 0. ENDIF zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN IF (t(i,k) .GE. t_glace_min) THEN zalpha_tr = a_tr_sca(3) ELSE zalpha_tr = a_tr_sca(4) ENDIF zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i)) pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi) frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi ! Nucleation with a factor of -1 instead of -0.5 zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i)) pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi) ENDIF ENDDO ! Scavenging through impaction in the underlying layer DO kk = k-1, 1, -1 DO i = 1, klon IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN IF (t(i,kk) .GE. t_glace_min) THEN zalpha_tr = a_tr_sca(1) ELSE zalpha_tr = a_tr_sca(2) ENDIF zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i)) pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi) frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi ENDIF ENDDO ENDDO !--save some variables for ice supersaturation ! DO i = 1, klon ! for memory rneb_seri(i,k) = rneb(i,k) ! for diagnostics rnebclr(i,k) = 1.0 - rneb(i,k) - rnebss(i,k) qvc(i,k) = zqs(i) * rneb(i,k) qclr(i,k) = MAX(1.e-10,zq(i) - qvc(i,k) - qss(i,k)) !--added by OB because of pathologic cases with the lognormal qcld(i,k) = qvc(i,k) + zcond(i) ENDDO !q_sat CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,1,.false.,zqsatl(:,k),zdqs(:)) CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,2,.false.,zqsats(:,k),zdqs(:)) ENDDO !====================================================================== ! END OF VERTICAL LOOP !====================================================================== ! Rain or snow at the surface (depending on the first layer temperature) DO i = 1, klon snow(i) = zifl(i) rain(i) = zrfl(i) ENDDO IF (ncoreczq>0) THEN WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.' ENDIF END SUBROUTINE LSCP !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ END MODULE LSCP_MOD