Index: LMDZ6/trunk/DefLists/field_def_lmdz.xml =================================================================== --- LMDZ6/trunk/DefLists/field_def_lmdz.xml (revision 5005) +++ LMDZ6/trunk/DefLists/field_def_lmdz.xml (revision 5007) @@ -634,4 +634,6 @@ + + @@ -642,5 +644,20 @@ - + + + + + + + + + + + + + + + + @@ -687,5 +704,8 @@ - + + + + Index: LMDZ6/trunk/libf/phylmd/lmdz_lscp.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lmdz_lscp.F90 (revision 5005) +++ LMDZ6/trunk/libf/phylmd/lmdz_lscp.F90 (revision 5007) @@ -7,13 +7,14 @@ !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ SUBROUTINE lscp(klon,klev,dtime,missing_val, & - paprs,pplay,temp,qt,ptconv,ratqs, & + paprs,pplay,temp,qt,qice_save,ptconv,ratqs, & d_t, d_q, d_ql, d_qi, rneb, rneblsvol, rneb_seri, & - pfraclr,pfracld, & + pfraclr, pfracld, & + cldfraliq, sigma2_icefracturb,mean_icefracturb, & radocond, radicefrac, rain, snow, & frac_impa, frac_nucl, beta, & - prfl, psfl, rhcl, qta, fraca, & - tv, pspsk, tla, thl, iflag_cld_th, & + prfl, psfl, rhcl, qta, fraca, & + tv, pspsk, tla, thl, iflag_cld_th, & iflag_ice_thermo, ok_ice_sursat, qsatl, qsats, & - distcltop,temp_cltop, & + distcltop, temp_cltop, tke, tke_dissip, & qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, & Tcontr, qcontr, qcontr2, fcontrN, fcontrP, & @@ -96,5 +97,6 @@ ! USE de modules contenant des fonctions. USE lmdz_cloudth, ONLY : cloudth, cloudth_v3, cloudth_v6, cloudth_mpc -USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, icefrac_lscp, calc_gammasat +USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, calc_gammasat +USE lmdz_lscp_tools, ONLY : icefrac_lscp, icefrac_lscp_turb USE lmdz_lscp_tools, ONLY : fallice_velocity, distance_to_cloud_top USE ice_sursat_mod, ONLY : ice_sursat @@ -111,5 +113,5 @@ USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG USE lmdz_lscp_ini, ONLY : ok_poprecip - +USE lmdz_lscp_ini, ONLY : iflag_icefrac IMPLICIT NONE @@ -129,4 +131,5 @@ REAL, DIMENSION(klon,klev), INTENT(IN) :: temp ! temperature (K) REAL, DIMENSION(klon,klev), INTENT(IN) :: qt ! total specific humidity (in vapor phase in input) [kg/kg] + REAL, DIMENSION(klon,klev), INTENT(IN) :: qice_save ! ice specific from previous time step [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 @@ -138,67 +141,72 @@ !Inputs associated with thermal plumes - REAL, DIMENSION(klon,klev), INTENT(IN) :: tv ! virtual potential temperature [K] - REAL, DIMENSION(klon,klev), INTENT(IN) :: qta ! 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) :: pspsk ! exner potential (p/100000)**(R/cp) - REAL, DIMENSION(klon,klev), INTENT(IN) :: tla ! liquid temperature within thermals [K] + REAL, DIMENSION(klon,klev), INTENT(IN) :: tv ! virtual potential temperature [K] + REAL, DIMENSION(klon,klev), INTENT(IN) :: qta ! 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) :: pspsk ! exner potential (p/100000)**(R/cp) + REAL, DIMENSION(klon,klev), INTENT(IN) :: tla ! liquid temperature within thermals [K] + REAL, DIMENSION(klon,klev), INTENT(IN) :: tke !--turbulent kinetic energy [m2/s2] + REAL, DIMENSION(klon,klev), INTENT(IN) :: tke_dissip !--TKE dissipation [m2/s3] ! INPUT/OUTPUT variables !------------------------ - REAL, DIMENSION(klon,klev), INTENT(INOUT) :: thl ! liquid potential temperature [K] - REAL, DIMENSION(klon,klev), INTENT(INOUT) :: ratqs ! function of pressure that sets the large-scale + REAL, DIMENSION(klon,klev), INTENT(INOUT) :: thl ! liquid potential temperature [K] + 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 + 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) :: rneblsvol ! cloud fraction per unit volume [-] - REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfraclr ! precip fraction clear-sky part [-] - REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfracld ! precip fraction cloudy part [-] - REAL, DIMENSION(klon,klev), INTENT(OUT) :: radocond ! 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 ! surface large-scale rainfall [kg/s/m2] - REAL, DIMENSION(klon), INTENT(OUT) :: snow ! surface large-scale snowfall [kg/s/m2] - REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsatl ! saturation specific humidity wrt liquid [kg/kg] - REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsats ! saturation specific humidity wrt ice [kg/kg] - 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] - REAL, DIMENSION(klon,klev), INTENT(OUT) :: distcltop ! distance to cloud top [m] - REAL, DIMENSION(klon,klev), INTENT(OUT) :: temp_cltop ! temperature of cloud top [K] - REAL, DIMENSION(klon,klev), INTENT(OUT) :: beta ! conversion rate of condensed water - - ! fraction of aerosol scavenging through impaction and nucleation (for on-line) + 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) :: rneblsvol ! cloud fraction per unit volume [-] + REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfraclr ! precip fraction clear-sky part [-] + REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfracld ! precip fraction cloudy part [-] + REAL, DIMENSION(klon,klev), INTENT(OUT) :: cldfraliq ! liquid fraction of cloud [-] + REAL, DIMENSION(klon,klev), INTENT(OUT) :: sigma2_icefracturb ! Variance of the diagnostic supersaturation distribution (icefrac_turb) [-] + REAL, DIMENSION(klon,klev), INTENT(OUT) :: mean_icefracturb ! Mean of the diagnostic supersaturation distribution (icefrac_turb) [-] + REAL, DIMENSION(klon,klev), INTENT(OUT) :: radocond ! 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 ! surface large-scale rainfall [kg/s/m2] + REAL, DIMENSION(klon), INTENT(OUT) :: snow ! surface large-scale snowfall [kg/s/m2] + REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsatl ! saturation specific humidity wrt liquid [kg/kg] + REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsats ! saturation specific humidity wrt ice [kg/kg] + 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] + REAL, DIMENSION(klon,klev), INTENT(OUT) :: distcltop ! distance to cloud top [m] + REAL, DIMENSION(klon,klev), INTENT(OUT) :: temp_cltop ! temperature of cloud top [K] + REAL, DIMENSION(klon,klev), INTENT(OUT) :: beta ! conversion rate of condensed water + + ! 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 [-] + 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 [-] ! for supersaturation and contrails parameterisation - REAL, DIMENSION(klon,klev), INTENT(OUT) :: qclr ! specific total water content in clear sky region [kg/kg] - REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcld ! specific total water content in cloudy region [kg/kg] - REAL, DIMENSION(klon,klev), INTENT(OUT) :: qss ! specific total water content in supersat region [kg/kg] - REAL, DIMENSION(klon,klev), INTENT(OUT) :: qvc ! specific vapor content in clouds [kg/kg] - REAL, DIMENSION(klon,klev), INTENT(OUT) :: rnebclr ! mesh fraction of clear sky [-] - REAL, DIMENSION(klon,klev), INTENT(OUT) :: rnebss ! mesh fraction of ISSR [-] - REAL, DIMENSION(klon,klev), INTENT(OUT) :: gamma_ss ! coefficient governing the ice nucleation RHi threshold [-] - REAL, DIMENSION(klon,klev), INTENT(OUT) :: Tcontr ! threshold temperature for contrail formation [K] - REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcontr ! threshold humidity for contrail formation [kg/kg] - REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcontr2 ! // (2nd expression more consistent with LMDZ expression of q) - REAL, DIMENSION(klon,klev), INTENT(OUT) :: fcontrN ! fraction of grid favourable to non-persistent contrails - REAL, DIMENSION(klon,klev), INTENT(OUT) :: fcontrP ! fraction of grid favourable to persistent contrails - REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sth ! mean saturation deficit in thermals - REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_senv ! mean saturation deficit in environment - REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sigmath ! std of saturation deficit in thermals - REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sigmaenv ! std of saturation deficit in environment + REAL, DIMENSION(klon,klev), INTENT(OUT) :: qclr ! specific total water content in clear sky region [kg/kg] + REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcld ! specific total water content in cloudy region [kg/kg] + REAL, DIMENSION(klon,klev), INTENT(OUT) :: qss ! specific total water content in supersat region [kg/kg] + REAL, DIMENSION(klon,klev), INTENT(OUT) :: qvc ! specific vapor content in clouds [kg/kg] + REAL, DIMENSION(klon,klev), INTENT(OUT) :: rnebclr ! mesh fraction of clear sky [-] + REAL, DIMENSION(klon,klev), INTENT(OUT) :: rnebss ! mesh fraction of ISSR [-] + REAL, DIMENSION(klon,klev), INTENT(OUT) :: gamma_ss ! coefficient governing the ice nucleation RHi threshold [-] + REAL, DIMENSION(klon,klev), INTENT(OUT) :: Tcontr ! threshold temperature for contrail formation [K] + REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcontr ! threshold humidity for contrail formation [kg/kg] + REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcontr2 ! // (2nd expression more consistent with LMDZ expression of q) + REAL, DIMENSION(klon,klev), INTENT(OUT) :: fcontrN ! fraction of grid favourable to non-persistent contrails + REAL, DIMENSION(klon,klev), INTENT(OUT) :: fcontrP ! fraction of grid favourable to persistent contrails + REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sth ! mean saturation deficit in thermals + REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_senv ! mean saturation deficit in environment + REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sigmath ! std of saturation deficit in thermals + REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sigmaenv ! std of saturation deficit in environment ! for POPRECIP @@ -221,6 +229,5 @@ ! LOCAL VARIABLES: !---------------- - - REAL,DIMENSION(klon) :: qsl, qsi + REAL,DIMENSION(klon) :: qsl, qsi ! saturation threshold at current vertical level REAL :: zct, zcl,zexpo REAL, DIMENSION(klon,klev) :: ctot @@ -229,10 +236,10 @@ REAL :: zdelta, zcor, zcvm5 REAL, DIMENSION(klon) :: zdqsdT_raw - REAL, DIMENSION(klon) :: gammasat,dgammasatdt ! coefficient to make cold condensation at the correct RH and derivative wrt T - REAL, DIMENSION(klon) :: Tbef,qlbef,DT + REAL, DIMENSION(klon) :: gammasat,dgammasatdt ! coefficient to make cold condensation at the correct RH and derivative wrt T + REAL, DIMENSION(klon) :: Tbef,qlbef,DT ! temperature, humidity and temp. variation during lognormal iteration REAL :: num,denom REAL :: cste - REAL, DIMENSION(klon) :: zpdf_sig,zpdf_k,zpdf_delta - REAL, DIMENSION(klon) :: Zpdf_a,zpdf_b,zpdf_e1,zpdf_e2 + REAL, DIMENSION(klon) :: zpdf_sig,zpdf_k,zpdf_delta ! lognormal parameters + REAL, DIMENSION(klon) :: Zpdf_a,zpdf_b,zpdf_e1,zpdf_e2 ! lognormal intermediate variables REAL :: erf REAL, DIMENSION(klon) :: zfice_th @@ -251,4 +258,5 @@ REAL :: zmelt,zrain,zsnow,zprecip REAL, DIMENSION(klon) :: dzfice + REAL, DIMENSION(klon) :: zfice_turb, dzfice_turb REAL :: zsolid REAL, DIMENSION(klon) :: qtot, qzero @@ -281,6 +289,6 @@ REAL, DIMENSION(klon,klev) :: radocondi, radocondl REAL :: effective_zneb - REAL, DIMENSION(klon) :: distcltop1D, temp_cltop1D - + REAL, DIMENSION(klon) :: zdistcltop, ztemp_cltop + REAL, DIMENSION(klon) :: zqliq, zqice, zqvapcl ! for icefrac_lscp_turb INTEGER i, k, n, kk, iter @@ -337,4 +345,7 @@ pfraclr(:,:)=0.0 pfracld(:,:)=0.0 +cldfraliq(:,:)=0. +sigma2_icefracturb(:,:)=0. +mean_icefracturb(:,:)=0. radocond(:,:) = 0.0 radicefrac(:,:) = 0.0 @@ -346,4 +357,6 @@ zfice(:)=0.0 dzfice(:)=0.0 +zfice_turb(:)=0.0 +dzfice_turb(:)=0.0 zqprecl(:)=0.0 zqpreci(:)=0.0 @@ -360,6 +373,6 @@ d_tot_zneb(:) = 0.0 qzero(:) = 0.0 -distcltop1D(:)=0.0 -temp_cltop1D(:) = 0.0 +zdistcltop(:)=0.0 +ztemp_cltop(:) = 0.0 ztupnew(:)=0.0 !--ice supersaturation @@ -374,4 +387,5 @@ distcltop(:,:)=0. temp_cltop(:,:)=0. + !-- poprecip qraindiag(:,:)= 0. @@ -393,5 +407,4 @@ - !c_iso: variable initialisation for iso @@ -412,4 +425,9 @@ ! Initialisation temperature and specific humidity + ! temp(klon,klev) is not modified by the routine, instead all changes in temperature are made on zt + ! at the end of the klon loop, a temperature incremtent d_t due to all processes + ! (thermalization, evap/sub incoming precip, cloud formation, precipitation processes) is calculated + ! d_t = temperature tendency due to lscp + ! The temperature of the overlying layer is updated here because needed for thermalization DO i = 1, klon zt(i)=temp(i,k) @@ -746,5 +764,5 @@ ELSEIF (iflag_cloudth_vert .EQ. 7) THEN ! Updated version of Arnaud Jam (correction by E. Vignon) + adapted treatment - ! for boundary-layer mixed phase clouds following Vignon et al. + ! for boundary-layer mixed phase clouds CALL cloudth_mpc(klon,klev,k,mpc_bl_points,zt,zq,qta(:,k),fraca(:,k), & pspsk(:,k),paprs(:,k+1),paprs(:,k),pplay(:,k), tla(:,k), & @@ -766,15 +784,15 @@ ! lognormal - lognormale = .TRUE. + lognormale(:) = .TRUE. ELSEIF (iflag_cld_th .GE. 6) THEN ! lognormal distribution when no thermals - lognormale = fraca(:,k) < min_frac_th_cld + lognormale(:) = fraca(:,k) < min_frac_th_cld ELSE ! When iflag_cld_th=5, we always assume ! bi-gaussian distribution - lognormale = .FALSE. + lognormale(:) = .FALSE. ENDIF @@ -829,7 +847,8 @@ IF (iflag_t_glace.GE.4) THEN ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top - CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop1D) + CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop) ENDIF - CALL icefrac_lscp(klon, zt(:), iflag_ice_thermo, distcltop1D(:),temp_cltop1D(:),zfice(:),dzfice(:)) + + CALL icefrac_lscp(klon, zt(:), iflag_ice_thermo, zdistcltop(:),ztemp_cltop(:),zfice(:),dzfice(:)) DO i=1,klon !todoan : check if loop in i is needed @@ -896,5 +915,6 @@ cste=RLSTT ENDIF - + + ! LEA_R : check formule qlbef(i)=max(0.,zqn(i)-zqs(i)) num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT & @@ -927,17 +947,25 @@ ! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top IF (iflag_t_glace.GE.4) THEN - CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,distcltop1D,temp_cltop1D) - distcltop(:,k)=distcltop1D(:) - temp_cltop(:,k)=temp_cltop1D(:) - ENDIF - ! Partition function in stratiform clouds (will be overwritten in boundary-layer MPCs) - CALL icefrac_lscp(klon,zt,iflag_ice_thermo,distcltop1D,temp_cltop1D,zfice,dzfice) - + CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop) + distcltop(:,k)=zdistcltop(:) + temp_cltop(:,k)=ztemp_cltop(:) + ENDIF + + ! Partition function depending on temperature + CALL icefrac_lscp(klon, Tbef, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice) + + ! Partition function depending on tke for non shallow-convective clouds + IF (iflag_icefrac .GE. 1) THEN + + CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), qice_save(:,k), ziflcld, zqn, & + rneb(:,k), tke(:,k), tke_dissip(:,k), zqliq, zqvapcl, zqice, zfice_turb, dzfice_turb, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k)) + + ENDIF ! Water vapor update, Phase determination and subsequent latent heat exchange DO i=1, klon - + ! Overwrite phase partitioning in boundary layer mixed phase clouds when the + ! iflag_cloudth_vert=7 and specific param is activated 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 @@ -946,7 +974,5 @@ zq(i) = zq(i) - zcond(i) zfice(i)=zfice_th(i) - ELSE - ! Checks on rneb, rhcl and zqn IF (rneb(i,k) .LE. 0.0) THEN @@ -964,4 +990,12 @@ ! following line is very strange and probably wrong: rhcl(i,k)=(zqs(i)+zq(i))/2./zqs(i) + ! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param) + IF (iflag_icefrac .GE. 1) THEN + IF (lognormale(i)) THEN + zcond(i) = zqliq(i) + zqice(i) + zfice(i)=zfice_turb(i) + rhcl(i,k) = zqvapcl(i) * rneb(i,k) + (zq(i) - zqn(i)) * (1.-rneb(i,k)) + ENDIF + ENDIF ENDIF @@ -1306,5 +1340,8 @@ znebprecipcld(i)=0.0 ENDIF - + !IF ( ((1-zfice(i))*zoliq(i) .GT. 0.) .AND. (zt(i) .LE. 233.15) ) THEN + !print*,'WARNING LEA OLIQ A <-40°C ' + !print*,'zt,Tbef,oliq,oice,cldfraliq,icefrac,rneb',zt(i),Tbef(i),(1-zfice(i))*zoliq(i),zfice(i)*zoliq(i),cldfraliq(i,k),zfice(i),rneb(i,k) + !ENDIF ENDDO Index: LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.F90 (revision 5005) +++ LMDZ6/trunk/libf/phylmd/lmdz_lscp_ini.F90 (revision 5007) @@ -6,6 +6,6 @@ !-------------------- - REAL RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RPI - !$OMP THREADPRIVATE(RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RPI) + REAL RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI + !$OMP THREADPRIVATE(RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI) REAL, SAVE, PROTECTED :: seuil_neb=0.001 ! cloud fraction threshold: a cloud can precipitate when exceeded @@ -67,29 +67,33 @@ !$OMP THREADPRIVATE(iflag_t_glace) - INTEGER, SAVE, PROTECTED :: iflag_cloudth_vert=0 ! option for determining cloud fraction and content in convective boundary layers + INTEGER, SAVE, PROTECTED :: iflag_cloudth_vert=0 ! option for determining cloud fraction and content in convective boundary layers !$OMP THREADPRIVATE(iflag_cloudth_vert) - INTEGER, SAVE, PROTECTED :: iflag_gammasat=0 ! which threshold for homogeneous nucleation below -40oC + INTEGER, SAVE, PROTECTED :: iflag_gammasat=0 ! which threshold for homogeneous nucleation below -40oC !$OMP THREADPRIVATE(iflag_gammasat) - INTEGER, SAVE, PROTECTED :: iflag_rain_incloud_vol=0 ! use of volume cloud fraction for rain autoconversion + INTEGER, SAVE, PROTECTED :: iflag_rain_incloud_vol=0 ! use of volume cloud fraction for rain autoconversion !$OMP THREADPRIVATE(iflag_rain_incloud_vol) - INTEGER, SAVE, PROTECTED :: iflag_bergeron=0 ! bergeron effect for liquid precipitation treatment + INTEGER, SAVE, PROTECTED :: iflag_bergeron=0 ! bergeron effect for liquid precipitation treatment !$OMP THREADPRIVATE(iflag_bergeron) - INTEGER, SAVE, PROTECTED :: iflag_fisrtilp_qsat=0 ! qsat adjustment (iterative) during autoconversion + INTEGER, SAVE, PROTECTED :: iflag_fisrtilp_qsat=0 ! qsat adjustment (iterative) during autoconversion !$OMP THREADPRIVATE(iflag_fisrtilp_qsat) - INTEGER, SAVE, PROTECTED :: iflag_pdf=0 ! type of subgrid scale qtot pdf + INTEGER, SAVE, PROTECTED :: iflag_pdf=0 ! type of subgrid scale qtot pdf !$OMP THREADPRIVATE(iflag_pdf) - INTEGER, SAVE, PROTECTED :: iflag_autoconversion=0 ! autoconversion option + INTEGER, SAVE, PROTECTED :: iflag_icefrac=0 ! which phase partitioning function to use + !$OMP THREADPRIVATE(iflag_icefrac) + + INTEGER, SAVE, PROTECTED :: iflag_autoconversion=0 ! autoconversion option !$OMP THREADPRIVATE(iflag_autoconversion) - LOGICAL, SAVE, PROTECTED :: reevap_ice=.false. ! no liquid precip for T< threshold + + LOGICAL, SAVE, PROTECTED :: reevap_ice=.false. ! no liquid precip for T< threshold !$OMP THREADPRIVATE(reevap_ice) - REAL, SAVE, PROTECTED :: cld_lc_lsc=2.6e-4 ! liquid autoconversion coefficient, stratiform rain + REAL, SAVE, PROTECTED :: cld_lc_lsc=2.6e-4 ! liquid autoconversion coefficient, stratiform rain !$OMP THREADPRIVATE(cld_lc_lsc) @@ -118,5 +122,5 @@ !$OMP THREADPRIVATE(coef_eva) - REAL, SAVE, PROTECTED :: coef_sub ! tuning coefficient ice precip sublimation + REAL, SAVE, PROTECTED :: coef_sub ! tuning coefficient ice precip sublimation !$OMP THREADPRIVATE(coef_sub) @@ -124,5 +128,5 @@ !$OMP THREADPRIVATE(expo_eva) - REAL, SAVE, PROTECTED :: expo_sub ! tuning coefficient ice precip sublimation + REAL, SAVE, PROTECTED :: expo_sub ! tuning coefficient ice precip sublimation !$OMP THREADPRIVATE(expo_sub) @@ -158,4 +162,22 @@ !$OMP THREADPRIVATE(thresh_precip_frac) + REAL, SAVE, PROTECTED :: tau_mixenv=100000 ! Homogeneization time of mixed phase clouds [s] + !$OMP THREADPRIVATE(tau_mixenv) + + REAL, SAVE, PROTECTED :: capa_crystal=1. ! Sursaturation of ice part in mixed phase clouds [-] + !$OMP THREADPRIVATE(capa_crystal) + + REAL, SAVE, PROTECTED :: lmix_mpc=1000 ! Length of turbulent zones in Mixed Phase Clouds [m] + !$OMP THREADPRIVATE(lmix_mpc) + + REAL, SAVE, PROTECTED :: naero5=0.5 ! Number concentration of aerosol larger than 0.5 microns [scm-3] + !$OMP THREADPRIVATE(naero5) + + REAL, SAVE, PROTECTED :: gamma_snwretro = 0. ! Proportion of snow taken into account in ice retroaction in icefrac_turb [-] + !$OMP THREADPRIVATE(gamma_snwretro) + + REAL, SAVE, PROTECTED :: gamma_taud = 1. ! Tuning coeff for tau_dissipturb [-] + !$OMP THREADPRIVATE(gamma_taud) + REAL, SAVE, PROTECTED :: gamma_col=1. ! A COMMENTER TODO [-] !$OMP THREADPRIVATE(gamma_col) @@ -167,17 +189,17 @@ !$OMP THREADPRIVATE(gamma_rim) - REAL, SAVE, PROTECTED :: rho_rain=1000. ! A COMMENTER TODO [kg/m3] + REAL, SAVE, PROTECTED :: rho_rain=1000. ! Rain density [kg/m3] !$OMP THREADPRIVATE(rho_rain) - REAL, SAVE, PROTECTED :: rho_ice=920. ! A COMMENTER TODO [kg/m3] + REAL, SAVE, PROTECTED :: rho_ice=920. ! Ice density [kg/m3] !$OMP THREADPRIVATE(rho_ice) - REAL, SAVE, PROTECTED :: r_rain=500.E-6 ! A COMMENTER TODO [m] + REAL, SAVE, PROTECTED :: r_rain=500.E-6 ! Rain droplets radius for POPRECIP [m] !$OMP THREADPRIVATE(r_rain) - REAL, SAVE, PROTECTED :: r_snow=1.E-3 ! A COMMENTER TODO [m] + REAL, SAVE, PROTECTED :: r_snow=1.E-3 ! Ice crystals radius for POPRECIP [m] !$OMP THREADPRIVATE(r_snow) - REAL, SAVE, PROTECTED :: tau_auto_snow_min=100. ! A COMMENTER TODO [s] + REAL, SAVE, PROTECTED :: tau_auto_snow_min=100. ! A COMMENTER TODO [s] !$OMP THREADPRIVATE(tau_auto_snow_min) @@ -188,32 +210,32 @@ !$OMP THREADPRIVATE(eps) - REAL, SAVE, PROTECTED :: gamma_melt=1. ! A COMMENTER TODO [-] + REAL, SAVE, PROTECTED :: gamma_melt=1. ! A COMMENTER TODO [-] !$OMP THREADPRIVATE(gamma_melt) - REAL, SAVE, PROTECTED :: alpha_freez=4. ! A COMMENTER TODO [-] + REAL, SAVE, PROTECTED :: alpha_freez=4. ! A COMMENTER TODO [-] !$OMP THREADPRIVATE(alpha_freez) - REAL, SAVE, PROTECTED :: beta_freez=0.1 ! A COMMENTER TODO [m-3.s-1] + REAL, SAVE, PROTECTED :: beta_freez=0.1 ! A COMMENTER TODO [m-3.s-1] !$OMP THREADPRIVATE(beta_freez) - REAL, SAVE, PROTECTED :: gamma_freez=1. ! A COMMENTER TODO [-] + REAL, SAVE, PROTECTED :: gamma_freez=1. ! A COMMENTER TODO [-] !$OMP THREADPRIVATE(gamma_freez) - REAL, SAVE, PROTECTED :: rain_fallspeed=4. ! A COMMENTER TODO [m/s] + REAL, SAVE, PROTECTED :: rain_fallspeed=4. ! A COMMENTER TODO [m/s] !$OMP THREADPRIVATE(rain_fallspeed) - REAL, SAVE, PROTECTED :: rain_fallspeed_clr ! A COMMENTER TODO [m/s] + REAL, SAVE, PROTECTED :: rain_fallspeed_clr ! A COMMENTER TODO [m/s] !$OMP THREADPRIVATE(rain_fallspeed_clr) - REAL, SAVE, PROTECTED :: rain_fallspeed_cld ! A COMMENTER TODO [m/s] + REAL, SAVE, PROTECTED :: rain_fallspeed_cld ! A COMMENTER TODO [m/s] !$OMP THREADPRIVATE(rain_fallspeed_cld) - REAL, SAVE, PROTECTED :: snow_fallspeed=1. ! A COMMENTER TODO [m/s] + REAL, SAVE, PROTECTED :: snow_fallspeed=1. ! A COMMENTER TODO [m/s] !$OMP THREADPRIVATE(snow_fallspeed) - REAL, SAVE, PROTECTED :: snow_fallspeed_clr ! A COMMENTER TODO [m/s] + REAL, SAVE, PROTECTED :: snow_fallspeed_clr ! A COMMENTER TODO [m/s] !$OMP THREADPRIVATE(snow_fallspeed_clr) - REAL, SAVE, PROTECTED :: snow_fallspeed_cld ! A COMMENTER TODO [m/s] + REAL, SAVE, PROTECTED :: snow_fallspeed_cld ! A COMMENTER TODO [m/s] !$OMP THREADPRIVATE(snow_fallspeed_cld) !--End of the parameters for poprecip @@ -226,5 +248,5 @@ SUBROUTINE lscp_ini(dtime,lunout_in,prt_level_in,ok_ice_sursat, iflag_ratqs, fl_cor_ebil_in, RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in, & - RVTMP2_in, RTT_in,RD_in,RG_in,RPI_in) + RVTMP2_in, RTT_in,RD_in,RG_in,RV_in,RPI_in) @@ -238,5 +260,5 @@ REAL, INTENT(IN) :: RCPD_in, RLSTT_in, RLVTT_in, RLMLT_in - REAL, INTENT(IN) :: RVTMP2_in, RTT_in, RD_in, RG_in, RPI_in + REAL, INTENT(IN) :: RVTMP2_in, RTT_in, RD_in, RG_in, RV_in, RPI_in character (len=20) :: modname='lscp_ini_mod' character (len=80) :: abort_message @@ -254,5 +276,5 @@ RLMLT=RLMLT_in RTT=RTT_in - RG=RG_in + RV=RV_in RVTMP2=RVTMP2_in RPI=RPI_in @@ -275,4 +297,5 @@ CALL getin_p('iflag_fisrtilp_qsat',iflag_fisrtilp_qsat) CALL getin_p('iflag_pdf',iflag_pdf) + CALL getin_p('iflag_icefrac',iflag_icefrac) CALL getin_p('reevap_ice',reevap_ice) CALL getin_p('cld_lc_lsc',cld_lc_lsc) @@ -296,4 +319,10 @@ CALL getin_p('dist_liq',dist_liq) CALL getin_p('tresh_cl',tresh_cl) + CALL getin_p('tau_mixenv',tau_mixenv) + CALL getin_p('capa_crystal',capa_crystal) + CALL getin_p('lmix_mpc',lmix_mpc) + CALL getin_p('naero5',naero5) + CALL getin_p('gamma_snwretro',gamma_snwretro) + CALL getin_p('gamma_taud',gamma_taud) CALL getin_p('iflag_oldbug_fisrtilp',iflag_oldbug_fisrtilp) CALL getin_p('ok_poprecip',ok_poprecip) @@ -334,4 +363,5 @@ WRITE(lunout,*) 'lscp_ini, iflag_fisrtilp_qsat:', iflag_fisrtilp_qsat WRITE(lunout,*) 'lscp_ini, iflag_pdf', iflag_pdf + WRITE(lunout,*) 'lscp_ini, iflag_icefrac', iflag_icefrac WRITE(lunout,*) 'lscp_ini, reevap_ice', reevap_ice WRITE(lunout,*) 'lscp_ini, cld_lc_lsc', cld_lc_lsc @@ -352,4 +382,10 @@ WRITE(lunout,*) 'lscp_ini, dist_liq', dist_liq WRITE(lunout,*) 'lscp_ini, tresh_cl', tresh_cl + WRITE(lunout,*) 'lscp_ini, tau_mixenv', tau_mixenv + WRITE(lunout,*) 'lscp_ini, capa_crystal', capa_crystal + WRITE(lunout,*) 'lscp_ini, lmix_mpc', lmix_mpc + WRITE(lunout,*) 'lscp_ini, naero5', naero5 + WRITE(lunout,*) 'lscp_ini, gamma_snwretro', gamma_snwretro + WRITE(lunout,*) 'lscp_ini, gamma_taud', gamma_taud WRITE(lunout,*) 'lscp_ini, iflag_oldbug_fisrtilp', iflag_oldbug_fisrtilp WRITE(lunout,*) 'lscp_ini, fl_cor_ebil', fl_cor_ebil Index: LMDZ6/trunk/libf/phylmd/lmdz_lscp_tools.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/lmdz_lscp_tools.F90 (revision 5005) +++ LMDZ6/trunk/libf/phylmd/lmdz_lscp_tools.F90 (revision 5007) @@ -136,5 +136,5 @@ CHARACTER (len = 80) :: abort_message - IF ((iflag_t_glace.LT.2) .OR. (iflag_t_glace.GT.6)) THEN + IF ((iflag_t_glace.LT.2)) THEN !.OR. (iflag_t_glace.GT.6)) THEN abort_message = 'lscp cannot be used if iflag_t_glace<2 or >6' CALL abort_physic(modname,abort_message,1) @@ -194,5 +194,5 @@ ! with CMIP6 function of temperature at cloud top - IF (iflag_t_glace .EQ. 5) THEN + IF ((iflag_t_glace .EQ. 5) .OR. (iflag_t_glace .EQ. 7)) THEN liqfrac_tmp = (temp(i)-t_glace_min) / (t_glace_max-t_glace_min) liqfrac_tmp = MIN(MAX(liqfrac_tmp,0.0),1.0) @@ -232,8 +232,7 @@ ENDIF ENDIF - + ENDDO ! klon - RETURN @@ -241,4 +240,269 @@ !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +SUBROUTINE ICEFRAC_LSCP_TURB(klon, dtime, temp, pplay, paprsdn, paprsup, qice_ini, snowcld, qtot_incl, cldfra, tke, tke_dissip, qliq, qvap_cld, qice, icefrac, dicefracdT, cldfraliq, sigma2_icefracturb, mean_icefracturb) +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + ! Compute the liquid, ice and vapour content (+ice fraction) based + ! on turbulence (see Fields 2014, Furtado 2016) +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + + + USE lmdz_lscp_ini, ONLY : prt_level, lunout + USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG, RV, RPI + USE lmdz_lscp_ini, ONLY : seuil_neb, temp_nowater + USE lmdz_lscp_ini, ONLY : tau_mixenv, lmix_mpc, naero5, gamma_snwretro, gamma_taud, capa_crystal + USE lmdz_lscp_ini, ONLY : eps + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: klon !--number of horizontal grid points + REAL, INTENT(IN) :: dtime !--time step [s] + + REAL, INTENT(IN), DIMENSION(klon) :: temp !--temperature + REAL, INTENT(IN), DIMENSION(klon) :: pplay !--pressure in the middle of the layer [Pa] + REAL, INTENT(IN), DIMENSION(klon) :: paprsdn !--pressure at the bottom interface of the layer [Pa] + REAL, INTENT(IN), DIMENSION(klon) :: paprsup !--pressure at the top interface of the layer [Pa] + REAL, INTENT(IN), DIMENSION(klon) :: qtot_incl !--specific total cloud water content, in-cloud content [kg/kg] + REAL, INTENT(IN), DIMENSION(klon) :: cldfra !--cloud fraction in gridbox [-] + REAL, INTENT(IN), DIMENSION(klon) :: tke !--turbulent kinetic energy [m2/s2] + REAL, INTENT(IN), DIMENSION(klon) :: tke_dissip !--TKE dissipation [m2/s3] + + REAL, INTENT(IN), DIMENSION(klon) :: qice_ini !--initial specific ice content gridbox-mean [kg/kg] + REAL, INTENT(IN), DIMENSION(klon) :: snowcld + REAL, INTENT(OUT), DIMENSION(klon) :: qliq !--specific liquid content gridbox-mean [kg/kg] + REAL, INTENT(OUT), DIMENSION(klon) :: qvap_cld !--specific cloud vapor content, gridbox-mean [kg/kg] + REAL, INTENT(OUT), DIMENSION(klon) :: qice !--specific ice content gridbox-mean [kg/kg] + REAL, INTENT(OUT), DIMENSION(klon) :: icefrac !--fraction of ice in condensed water [-] + REAL, INTENT(OUT), DIMENSION(klon) :: dicefracdT + + REAL, INTENT(OUT), DIMENSION(klon) :: cldfraliq !--fraction of cldfra which is liquid only + REAL, INTENT(OUT), DIMENSION(klon) :: sigma2_icefracturb !--Temporary + REAL, INTENT(OUT), DIMENSION(klon) :: mean_icefracturb !--Temporary + + REAL, DIMENSION(klon) :: qzero, qsatl, dqsatl, qsati, dqsati !--specific humidity saturation values + INTEGER :: i + + REAL :: qvap_incl, qice_incl, qliq_incl, qiceini_incl !--In-cloud specific quantities [kg/kg] + REAL :: qsnowcld_incl + !REAL :: capa_crystal !--Capacitance of ice crystals [-] + REAL :: water_vapor_diff !--Water-vapour diffusion coefficient in air [m2/s] (function of T&P) + REAL :: air_thermal_conduct !--Thermal conductivity of air [J/m/K/s] (function of T) + REAL :: C0 !--Lagrangian structure function [-] + REAL :: tau_mixingenv + REAL :: tau_dissipturb + REAL :: invtau_phaserelax + REAL :: sigma2_pdf, mean_pdf + REAL :: ai, bi, B0 + REAL :: sursat_iceliq + REAL :: sursat_env + REAL :: liqfra_max + REAL :: sursat_iceext + REAL :: nb_crystals !--number concentration of ice crystals [#/m3] + REAL :: moment1_PSD !--1st moment of ice PSD + REAL :: N0_PSD, lambda_PSD !--parameters of the exponential PSD + + REAL :: rho_ice !--ice density [kg/m3] + REAL :: cldfra1D + REAL :: deltaz, rho_air + REAL :: psati !--saturation vapor pressure wrt i [Pa] + + C0 = 10. !--value assumed in Field2014 + rho_ice = 950. + sursat_iceext = -0.1 + !capa_crystal = 1. !r_ice + qzero(:) = 0. + cldfraliq(:) = 0. + icefrac(:) = 0. + dicefracdT(:) = 0. + + sigma2_icefracturb(:) = 0. + mean_icefracturb(:) = 0. + + !--wrt liquid water + CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,1,.false.,qsatl(:),dqsatl(:)) + !--wrt ice + CALL calc_qsat_ecmwf(klon,temp(:),qzero(:),pplay(:),RTT,2,.false.,qsati(:),dqsati(:)) + + + DO i=1,klon + + + rho_air = pplay(i) / temp(i) / RD + !deltaz = ( paprsdn(i) - paprsup(i) ) / RG / rho_air(i) + ! because cldfra is intent in, but can be locally modified due to test + cldfra1D = cldfra(i) + IF (cldfra(i) .LE. 0.) THEN + qvap_cld(i) = 0. + qliq(i) = 0. + qice(i) = 0. + cldfraliq(i) = 0. + icefrac(i) = 0. + dicefracdT(i) = 0. + + ! If there is a cloud + ELSE + IF (cldfra(i) .GE. 1.0) THEN + cldfra1D = 1.0 + END IF + + ! T>0°C, no ice allowed + IF ( temp(i) .GE. RTT ) THEN + qvap_cld(i) = qsatl(i) * cldfra1D + qliq(i) = MAX(0.0,qtot_incl(i)-qsatl(i)) * cldfra1D + qice(i) = 0. + cldfraliq(i) = 1. + icefrac(i) = 0. + dicefracdT(i) = 0. + + ! T<-38°C, no liquid allowed + ELSE IF ( temp(i) .LE. temp_nowater) THEN + qvap_cld(i) = qsati(i) * cldfra1D + qliq(i) = 0. + qice(i) = MAX(0.0,qtot_incl(i)-qsati(i)) * cldfra1D + cldfraliq(i) = 0. + icefrac(i) = 1. + dicefracdT(i) = 0. + + ! MPC temperature + ELSE + ! Not enough TKE + IF ( tke_dissip(i) .LE. eps ) THEN + qvap_cld(i) = qsati(i) * cldfra1D + qliq(i) = 0. + qice(i) = MAX(0.,qtot_incl(i)-qsati(i)) * cldfra1D + cldfraliq(i) = 0. + icefrac(i) = 1. + dicefracdT(i) = 0. + + ! Enough TKE + ELSE + !--------------------------------------------------------- + !-- ICE SUPERSATURATION PDF + !--------------------------------------------------------- + !--If -38°C< T <0°C and there is enough turbulence, + !--we compute the cloud liquid properties with a Gaussian PDF + !--of ice supersaturation F(Si) (Field2014, Furtado2016). + !--Parameters of the PDF are function of turbulence and + !--microphysics/existing ice. + + sursat_iceliq = qsatl(i)/qsati(i) - 1. + psati = qsati(i) * pplay(i) / (RD/RV) + + !-------------- MICROPHYSICAL TERMS -------------- + !--We assume an exponential ice PSD whose parameters + !--are computed following Morrison&Gettelman 2008 + !--Ice number density is assumed equals to INP density + !--which is a function of temperature (DeMott 2010) + !--bi and B0 are microphysical function characterizing + !--vapor/ice interactions + !--tau_phase_relax is the typical time of vapor deposition + !--onto ice crystals + + qiceini_incl = qice_ini(i) / cldfra1D + qsnowcld_incl = snowcld(i) * RG * dtime / ( paprsdn(i) - paprsup(i) ) / cldfra1D + sursat_env = max(0., (qtot_incl(i) - qiceini_incl)/qsati(i) - 1.) + IF ( qiceini_incl .GT. eps ) THEN + nb_crystals = 1.e3 * 5.94e-5 * ( RTT - temp(i) )**3.33 * naero5**(0.0264*(RTT-temp(i))+0.0033) + lambda_PSD = ( (RPI*rho_ice*nb_crystals*24.) / (6.*(qiceini_incl + gamma_snwretro * qsnowcld_incl)) ) ** (1./3.) + N0_PSD = nb_crystals * lambda_PSD + moment1_PSD = N0_PSD/2./lambda_PSD**2 + ELSE + moment1_PSD = 0. + ENDIF + + !--Formulae for air thermal conductivity and water vapor diffusivity + !--comes respectively from Beard and Pruppacher (1971) + !--and Hall and Pruppacher (1976) + + air_thermal_conduct = ( 5.69 + 0.017 * ( temp(i) - RTT ) ) * 1.e-3 * 4.184 + water_vapor_diff = 2.11*1e-5 * ( temp(i) / RTT )**1.94 * ( 101325 / pplay(i) ) + + bi = 1./((qsati(i)+qsatl(i))/2.) + RLSTT**2 / RCPD / RV / temp(i)**2 + B0 = 4. * RPI * capa_crystal * 1. / ( RLSTT**2 / air_thermal_conduct / RV / temp(i)**2 & + + RV * temp(i) / psati / water_vapor_diff ) + + invtau_phaserelax = (bi * B0 * moment1_PSD ) + +! Old way of estimating moment1 : spherical crystals + monodisperse PSD +! nb_crystals = rho_air * qiceini_incl / ( 4. / 3. * RPI * r_ice**3. * rho_ice ) +! moment1_PSD = nb_crystals * r_ice + + !----------------- TURBULENT SOURCE/SINK TERMS ----------------- + !--Tau_mixingenv is the time needed to homogeneize the parcel + !--with its environment by turbulent diffusion over the parcel + !--length scale + !--if lmix_mpc <0, tau_mixigenv value is prescribed + !--else tau_mixigenv value is derived from tke_dissip and lmix_mpc + !--Tau_dissipturb is the time needed turbulence to decay due to + !--viscosity + + ai = RG / RD / temp(i) * ( RD * RLSTT / RCPD / RV / temp(i) - 1. ) + IF ( lmix_mpc .GT. 0 ) THEN + tau_mixingenv = ( lmix_mpc**2. / tke_dissip(i) )**(1./3.) + ELSE + tau_mixingenv = tau_mixenv + ENDIF + + tau_dissipturb = gamma_taud * 2. * 2./3. * tke(i) / tke_dissip(i) / C0 + + !--------------------- PDF COMPUTATIONS --------------------- + !--Formulae for sigma2_pdf (variance), mean of PDF in Furtado2016 + !--cloud liquid fraction and in-cloud liquid content are given + !--by integrating resp. F(Si) and Si*F(Si) + !--Liquid is limited by the available water vapor trough a + !--maximal liquid fraction + + liqfra_max = MAX(0., (MIN (1.,( qtot_incl(i) - qiceini_incl - qsati(i) * (1 + sursat_iceext ) ) / ( qsatl(i) - qsati(i) ) ) ) ) + sigma2_pdf = 1./2. * ( ai**2 ) * 2./3. * tke(i) * tau_dissipturb / ( invtau_phaserelax + 1./tau_mixingenv ) + mean_pdf = sursat_env * 1./tau_mixingenv / ( invtau_phaserelax + 1./tau_mixingenv ) + cldfraliq(i) = 0.5 * (1. - erf( ( sursat_iceliq - mean_pdf) / (SQRT(2.* sigma2_pdf) ) ) ) + !IF (cldfraliq(i) .GT. liqfra_max) THEN + ! cldfraliq(i) = liqfra_max + !ENDIF + + qliq_incl = qsati(i) * SQRT(sigma2_pdf) / SQRT(2.*RPI) * EXP( -1.*(sursat_iceliq - mean_pdf)**2. / (2.*sigma2_pdf) ) & + - qsati(i) * cldfraliq(i) * (sursat_iceliq - mean_pdf ) + + sigma2_icefracturb(i)= sigma2_pdf + mean_icefracturb(i) = mean_pdf + !------------ ICE AMOUNT AND WATER CONSERVATION ------------ + + IF ( (qliq_incl .LE. eps) .OR. (cldfraliq(i) .LE. eps) ) THEN + qliq_incl = 0. + cldfraliq(i) = 0. + END IF + + !--Choice for in-cloud vapor : + !--1.Weighted mean between qvap in MPC parts and in ice-only parts + !--2.Always at ice saturation + qvap_incl = MAX(qsati(i), ( 1. - cldfraliq(i) ) * (sursat_iceext + 1.) * qsati(i) + cldfraliq(i) * qsatl(i) ) + + IF ( qvap_incl .GE. qtot_incl(i) ) THEN + qvap_incl = qsati(i) + qliq_incl = qtot_incl(i) - qvap_incl + qice_incl = 0. + + ELSEIF ( (qvap_incl + qliq_incl) .GE. qtot_incl(i) ) THEN + qliq_incl = MAX(0.0,qtot_incl(i) - qvap_incl) + qice_incl = 0. + ELSE + qice_incl = qtot_incl(i) - qvap_incl - qliq_incl + END IF + + qvap_cld(i) = qvap_incl * cldfra1D + qliq(i) = qliq_incl * cldfra1D + qice(i) = qice_incl * cldfra1D + icefrac(i) = qice(i) / ( qice(i) + qliq(i) ) + dicefracdT(i) = 0. + !print*,'MPC turb' + + END IF ! Enough TKE + + END IF ! ! MPC temperature + + END IF ! cldfra + + ENDDO ! klon +END SUBROUTINE ICEFRAC_LSCP_TURB +!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Index: LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90 (revision 5005) +++ LMDZ6/trunk/libf/phylmd/phys_local_var_mod.F90 (revision 5007) @@ -478,4 +478,10 @@ REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: pfraclr,pfracld !$OMP THREADPRIVATE(pfraclr,pfracld) + REAL, SAVE, ALLOCATABLE :: cldfraliq(:,:) +!$OMP THREADPRIVATE(cldfraliq) + REAL, SAVE, ALLOCATABLE ::mean_icefracturb(:,:) +!$OMP THREADPRIVATE(mean_icefracturb) + REAL, SAVE, ALLOCATABLE :: sigma2_icefracturb(:,:) +!$OMP THREADPRIVATE(sigma2_icefracturb) ! variables de sorties MM @@ -957,4 +963,7 @@ ALLOCATE(pfraclr(klon,klev),pfracld(klon,klev)) pfraclr(:,:)=0. ; pfracld(:,:)=0. ! because not always defined + ALLOCATE(cldfraliq(klon,klev)) + ALLOCATE(sigma2_icefracturb(klon,klev)) + ALLOCATE(mean_icefracturb(klon,klev)) ALLOCATE(distcltop(klon,klev)) ALLOCATE(temp_cltop(klon,klev)) @@ -1283,4 +1292,7 @@ DEALLOCATE(rneb) DEALLOCATE(pfraclr,pfracld) + DEALLOCATE(cldfraliq) + DEALLOCATE(sigma2_icefracturb) + DEALLOCATE(mean_icefracturb) DEALLOCATE (zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic) DEALLOCATE(distcltop) Index: LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90 (revision 5005) +++ LMDZ6/trunk/libf/phylmd/phys_output_ctrlout_mod.F90 (revision 5007) @@ -1559,4 +1559,11 @@ TYPE(ctrl_out), SAVE :: o_rneb = ctrl_out((/ 2, 5, 10, 10, 10, 10, 11, 11, 11, 11/), & 'rneb', 'Cloud fraction', '-', (/ ('', i=1, 10) /)) + TYPE(ctrl_out), SAVE :: o_cldfraliq = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), & + 'cldfraliq', 'Liquid fraction of the cloud', '-', (/ ('', i=1, 10) /)) + TYPE(ctrl_out), SAVE :: o_sigma2_icefracturb = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), & + 'sigma2_icefracturb', 'Variance of the diagnostic supersaturation distribution (icefrac_turb) [-]', '-', (/ ('', i=1, 10) /)) + TYPE(ctrl_out), SAVE :: o_mean_icefracturb = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), & + 'mean_icefracturb', 'Mean of the diagnostic supersaturation distribution (icefrac_turb) [-]', '-', (/ ('', i=1, 10) /)) + TYPE(ctrl_out), SAVE :: o_rnebjn = ctrl_out((/ 2, 5, 10, 10, 10, 10, 11, 11,11, 11/), & 'rnebjn', 'Cloud fraction in day', '-', (/ ('', i=1, 10) /)) Index: LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90 (revision 5005) +++ LMDZ6/trunk/libf/phylmd/phys_output_write_mod.F90 (revision 5007) @@ -140,5 +140,5 @@ o_zfull, o_zhalf, o_rneb, o_rnebjn, o_rnebcon, & o_rnebls, o_rneblsvol, o_rhum, o_rhl, o_rhi, o_ozone, o_ozone_light, & - o_pfraclr, o_pfracld, & + o_pfraclr, o_pfracld, o_cldfraliq, o_sigma2_icefracturb, o_mean_icefracturb, & o_qrainlsc, o_qsnowlsc, o_dqreva, o_dqrauto, o_dqrcol, o_dqrmelt, o_dqrfreez, & o_dqssub, o_dqsauto, o_dqsagg, o_dqsrim, o_dqsmelt, o_dqsfreez, & @@ -363,6 +363,7 @@ ql_seri, qs_seri, qbs_seri, tr_seri, qbs_seri,& zphi, u_seri, v_seri, omega, cldfra, & - rneb, rnebjn, rneblsvol, zx_rh, zx_rhl, zx_rhi, & - pfraclr, pfracld, & + rneb, rnebjn, rneblsvol, & + zx_rh, zx_rhl, zx_rhi, & + pfraclr, pfracld, cldfraliq, sigma2_icefracturb, mean_icefracturb, & qraindiag, qsnowdiag, dqreva, dqssub, & dqrauto,dqrcol,dqrmelt,dqrfreez, & @@ -1999,4 +2000,7 @@ CALL histwrite_phy(o_pfraclr, pfraclr) CALL histwrite_phy(o_pfracld, pfracld) + CALL histwrite_phy(o_cldfraliq, cldfraliq) + CALL histwrite_phy(o_sigma2_icefracturb, sigma2_icefracturb) + CALL histwrite_phy(o_mean_icefracturb, mean_icefracturb) IF (ok_poprecip) THEN CALL histwrite_phy(o_qrainlsc, qraindiag) Index: LMDZ6/trunk/libf/phylmd/physiq_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (revision 5005) +++ LMDZ6/trunk/libf/phylmd/physiq_mod.F90 (revision 5007) @@ -332,6 +332,6 @@ ! rneblsvol, & - pfraclr,pfracld, & - distcltop,temp_cltop, & + pfraclr, pfracld, cldfraliq, sigma2_icefracturb, mean_icefracturb, & + distcltop, temp_cltop, & zqsatl, zqsats, & qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, & @@ -1863,5 +1863,5 @@ & RG,RD,RCPD,RKAPPA,RLVTT,RETV) CALL ratqs_ini(klon,klev,iflag_thermals,lunout,nbsrf,is_lic,is_ter,RG,RV,RD,RCPD,RLSTT,RLVTT,RTT) - CALL lscp_ini(pdtphys,lunout,prt_level,ok_ice_sursat,iflag_ratqs,fl_cor_ebil,RCPD,RLSTT,RLVTT,RLMLT,RVTMP2,RTT,RD,RG,RPI) + CALL lscp_ini(pdtphys,lunout,prt_level,ok_ice_sursat,iflag_ratqs,fl_cor_ebil,RCPD,RLSTT,RLVTT,RLMLT,RVTMP2,RTT,RD,RG,RV,RPI) CALL blowing_snow_ini(RCPD, RLSTT, RLVTT, RLMLT, & RVTMP2, RTT,RD,RG, RV, RPI) @@ -1943,8 +1943,11 @@ ELSE IF (klon_glo==1) THEN pbl_tke(:,:,is_ave) = 0. + pbl_eps(:,:,is_ave) = 0. DO nsrf=1,nbsrf DO k = 1,klev+1 pbl_tke(:,k,is_ave) = pbl_tke(:,k,is_ave) & +pctsrf(:,nsrf)*pbl_tke(:,k,nsrf) + pbl_eps(:,k,is_ave) = pbl_eps(:,k,is_ave) & + +pctsrf(:,nsrf)*pbl_eps(:,k,nsrf) ENDDO ENDDO @@ -1952,4 +1955,5 @@ pbl_tke(:,:,is_ave) = 0. !ym missing init : maybe must be initialized in the same way that for klon_glo==1 ?? !>jyg + pbl_eps(:,:,is_ave) = 0. ENDIF !IM begin @@ -3896,13 +3900,13 @@ CALL lscp(klon,klev,phys_tstep,missing_val,paprs,pplay, & - t_seri, q_seri,ptconv,ratqs, & + t_seri, q_seri,qs_ancien,ptconv,ratqs, & d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, rneb_seri, & - pfraclr,pfracld, & + pfraclr, pfracld, cldfraliq, sigma2_icefracturb, mean_icefracturb, & radocond, picefra, rain_lsc, snow_lsc, & frac_impa, frac_nucl, beta_prec_fisrt, & prfl, psfl, rhcl, & zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, & - iflag_ice_thermo, ok_ice_sursat, zqsatl, zqsats, distcltop, temp_cltop, & - qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, & + iflag_ice_thermo, ok_ice_sursat, zqsatl, zqsats, distcltop, temp_cltop, & + pbl_tke(:,:,is_ave), pbl_eps(:,:,is_ave), qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, & Tcontr, qcontr, qcontr2, fcontrN, fcontrP , & cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, & Index: LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90 (revision 5005) +++ LMDZ6/trunk/libf/phylmdiso/phys_local_var_mod.F90 (revision 5007) @@ -608,4 +608,10 @@ REAL, ALLOCATABLE, SAVE, DIMENSION(:,:) :: pfraclr,pfracld !$OMP THREADPRIVATE(pfraclr,pfracld) + REAL, SAVE, ALLOCATABLE :: cldfraliq(:,:) +!$OMP THREADPRIVATE(cldfraliq) + REAL, SAVE, ALLOCATABLE ::mean_icefracturb(:,:) +!$OMP THREADPRIVATE(mean_icefracturb) + REAL, SAVE, ALLOCATABLE :: sigma2_icefracturb(:,:) +!$OMP THREADPRIVATE(sigma2_icefracturb) ! variables de sorties MM @@ -1156,4 +1162,7 @@ ALLOCATE(pfraclr(klon,klev),pfracld(klon,klev)) pfraclr(:,:)=0. ; pfracld(:,:)=0. ! because not always defined + ALLOCATE(cldfraliq(klon,klev)) + ALLOCATE(sigma2_icefracturb(klon,klev)) + ALLOCATE(mean_icefracturb(klon,klev)) ALLOCATE(distcltop(klon,klev)) ALLOCATE(temp_cltop(klon,klev)) @@ -1543,4 +1552,7 @@ DEALLOCATE(rneb) DEALLOCATE(pfraclr,pfracld) + DEALLOCATE(cldfraliq) + DEALLOCATE(sigma2_icefracturb) + DEALLOCATE(mean_icefracturb) DEALLOCATE (zxsnow,snowhgt,qsnow,to_ice,sissnow,runoff,albsol3_lic) DEALLOCATE(distcltop) Index: LMDZ6/trunk/libf/phylmdiso/phys_output_ctrlout_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/phys_output_ctrlout_mod.F90 (revision 5005) +++ LMDZ6/trunk/libf/phylmdiso/phys_output_ctrlout_mod.F90 (revision 5007) @@ -1553,4 +1553,11 @@ TYPE(ctrl_out), SAVE :: o_rneb = ctrl_out((/ 2, 5, 10, 10, 10, 10, 11, 11, 11, 11/), & 'rneb', 'Cloud fraction', '-', (/ ('', i=1, 10) /)) + TYPE(ctrl_out), SAVE :: o_cldfraliq = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), & + 'cldfraliq', 'Liquid fraction of the cloud', '-', (/ ('', i=1, 10) /)) + TYPE(ctrl_out), SAVE :: o_sigma2_icefracturb = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), & + 'sigma2_icefracturb', 'Variance of the diagnostic supersaturation distribution (icefrac_turb) [-]', '-', (/ ('', i=1, 10) /)) + TYPE(ctrl_out), SAVE :: o_mean_icefracturb = ctrl_out((/ 11, 11, 11, 11, 11, 11, 11, 11, 11, 11/), & + 'mean_icefracturb', 'Mean of the diagnostic supersaturation distribution (icefrac_turb) [-]', '-', (/ ('', i=1, 10) /)) + TYPE(ctrl_out), SAVE :: o_rnebjn = ctrl_out((/ 2, 5, 10, 10, 10, 10, 11, 11,11, 11/), & 'rnebjn', 'Cloud fraction in day', '-', (/ ('', i=1, 10) /)) Index: LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90 =================================================================== --- LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90 (revision 5005) +++ LMDZ6/trunk/libf/phylmdiso/physiq_mod.F90 (revision 5007) @@ -373,6 +373,6 @@ ! rneblsvol, & - pfraclr,pfracld, & - distcltop,temp_cltop, & + pfraclr, pfracld, cldfraliq, sigma2_icefracturb, mean_icefracturb, & + distcltop, temp_cltop, & zqsatl, zqsats, & qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, & @@ -2026,5 +2026,5 @@ & RG,RD,RCPD,RKAPPA,RLVTT,RETV) CALL ratqs_ini(klon,klev,iflag_thermals,lunout,nbsrf,is_lic,is_ter,RG,RV,RD,RCPD,RLSTT,RLVTT,RTT) - CALL lscp_ini(pdtphys,lunout,prt_level,ok_ice_sursat,iflag_ratqs,fl_cor_ebil,RCPD,RLSTT,RLVTT,RLMLT,RVTMP2,RTT,RD,RG,RPI) + CALL lscp_ini(pdtphys,lunout,prt_level,ok_ice_sursat,iflag_ratqs,fl_cor_ebil,RCPD,RLSTT,RLVTT,RLMLT,RVTMP2,RTT,RD,RG,RV,RPI) CALL blowing_snow_ini(RCPD, RLSTT, RLVTT, RLMLT, & RVTMP2, RTT,RD,RG, RV, RPI) @@ -5076,13 +5076,13 @@ CALL lscp(klon,klev,phys_tstep,missing_val,paprs,pplay, & - t_seri, q_seri,ptconv,ratqs, & + t_seri, q_seri,qs_ancien,ptconv,ratqs, & d_t_lsc, d_q_lsc, d_ql_lsc, d_qi_lsc, rneb, rneblsvol, rneb_seri, & - pfraclr,pfracld, & + pfraclr, pfracld, cldfraliq, sigma2_icefracturb, mean_icefracturb, & radocond, picefra, rain_lsc, snow_lsc, & frac_impa, frac_nucl, beta_prec_fisrt, & prfl, psfl, rhcl, & zqasc, fraca,ztv,zpspsk,ztla,zthl,iflag_cld_th, & - iflag_ice_thermo, ok_ice_sursat, zqsatl, zqsats, distcltop, temp_cltop, & - qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, & + iflag_ice_thermo, ok_ice_sursat, zqsatl, zqsats, distcltop, temp_cltop, & + pbl_tke(:,:,is_ave), pbl_eps(:,:,is_ave), qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, & Tcontr, qcontr, qcontr2, fcontrN, fcontrP , & cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &