Context Navigation

lmdz_lscp.F90 @ 5082

Last change on this file since 5082 was 5082, checked in by abarral, 4 months ago
(lint) Fix obsolete boolean operators
File size: 64.5 KB

Rev	Line
[4664]	1	MODULE lmdz_lscp
[3999]	2
	3	IMPLICIT NONE
	4
	5	CONTAINS
	6
	7	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[4380]	8	SUBROUTINE lscp(klon,klev,dtime,missing_val, &
[5007]	9	paprs,pplay,temp,qt,qice_save,ptconv,ratqs, &
[4380]	10	d_t, d_q, d_ql, d_qi, rneb, rneblsvol, rneb_seri, &
[5007]	11	pfraclr, pfracld, &
	12	cldfraliq, sigma2_icefracturb,mean_icefracturb, &
[4530]	13	radocond, radicefrac, rain, snow, &
[3999]	14	frac_impa, frac_nucl, beta, &
[5007]	15	prfl, psfl, rhcl, qta, fraca, &
	16	tv, pspsk, tla, thl, iflag_cld_th, &
[4380]	17	iflag_ice_thermo, ok_ice_sursat, qsatl, qsats, &
[5007]	18	distcltop, temp_cltop, tke, tke_dissip, &
[4380]	19	qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, &
[4651]	20	Tcontr, qcontr, qcontr2, fcontrN, fcontrP, &
[4803]	21	cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv, &
[4830]	22	qraindiag, qsnowdiag, dqreva, dqssub, dqrauto, &
	23	dqrcol, dqrmelt, dqrfreez, dqsauto, dqsagg, dqsrim,&
[4819]	24	dqsmelt, dqsfreez)
[3999]	25
	26	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	27	! Authors: Z.X. Li (LMD), J-L Dufresne (LMD), C. Rio (LMD), J-Y Grandpeix (LMD)
	28	! A. JAM (LMD), J-B Madeleine (LMD), E. Vignon (LMD), L. Touzze-Peiffert (LMD)
	29	!--------------------------------------------------------------------------------
	30	! Date: 01/2021
	31	!--------------------------------------------------------------------------------
	32	! Aim: Large Scale Clouds and Precipitation (LSCP)
[4226]	33	!
[3999]	34	! This code is a new version of the fisrtilp.F90 routine, which is itself a
[4072]	35	! merge of 'first' (superrsaturation physics, P. LeVan K. Laval)
[3999]	36	! and 'ilp' (il pleut, L. Li)
	37	!
[4380]	38	! Compared to the original fisrtilp code, lscp
[3999]	39	! -> assumes thermcep = .TRUE. all the time (fisrtilp inconsistent when .FALSE.)
	40	! -> consider always precipitation thermalisation (fl_cor_ebil>0)
	41	! -> option iflag_fisrtilp_qsat<0 no longer possible (qsat does not evolve with T)
	42	! -> option "oldbug" by JYG has been removed
	43	! -> iflag_t_glace >0 always
	44	! -> the 'all or nothing' cloud approach is no longer available (cpartiel=T always)
	45	! -> rectangular distribution from L. Li no longer available
	46	! -> We always account for the Wegener-Findeisen-Bergeron process (iflag_bergeron = 2 in fisrt)
	47	!--------------------------------------------------------------------------------
	48	! References:
	49	!
[4412]	50	! - Bony, S., & Emanuel, K. A. 2001, JAS, doi: 10.1175/1520-0469(2001)058<3158:APOTCA>2.0.CO;2
[3999]	51	! - Hourdin et al. 2013, Clim Dyn, doi:10.1007/s00382-012-1343-y
	52	! - Jam et al. 2013, Boundary-Layer Meteorol, doi:10.1007/s10546-012-9789-3
[4412]	53	! - Jouhaud, et al. 2018. JAMES, doi:10.1029/2018MS001379
[3999]	54	! - Madeleine et al. 2020, JAMES, doi:10.1029/2020MS002046
[4412]	55	! - Touzze-Peifert Ludo, PhD thesis, p117-124
[3999]	56	! -------------------------------------------------------------------------------
	57	! Code structure:
	58	!
[4226]	59	! P0> Thermalization of the precipitation coming from the overlying layer
[3999]	60	! P1> Evaporation of the precipitation (falling from the k+1 level)
	61	! P2> Cloud formation (at the k level)
[4412]	62	! P2.A.1> With the PDFs, calculation of cloud properties using the inital
[3999]	63	! values of T and Q
	64	! P2.A.2> Coupling between condensed water and temperature
	65	! P2.A.3> Calculation of final quantities associated with cloud formation
[4412]	66	! P2.B> Release of Latent heat after cloud formation
[3999]	67	! P3> Autoconversion to precipitation (k-level)
	68	! P4> Wet scavenging
	69	!------------------------------------------------------------------------------
	70	! Some preliminary comments (JBM) :
	71	!
	72	! The cloud water that the radiation scheme sees is not the same that the cloud
	73	! water used in the physics and the dynamics
	74	!
[4412]	75	! During the autoconversion to precipitation (P3 step), radocond (cloud water used
[3999]	76	! by the radiation scheme) is calculated as an average of the water that remains
	77	! in the cloud during the precipitation and not the water remaining at the end
	78	! of the time step. The latter is used in the rest of the physics and advected
	79	! by the dynamics.
	80	!
	81	! In summary:
	82	!
	83	! Radiation:
	84	! xflwc(newmicro)+xfiwc(newmicro) =
[4412]	85	! radocond=lwcon(nc)+iwcon(nc)
[3999]	86	!
	87	! Notetheless, be aware of:
	88	!
[4412]	89	! radocond .NE. ocond(nc)
[3999]	90	! i.e.:
	91	! lwcon(nc)+iwcon(nc) .NE. ocond(nc)
[4412]	92	! but oliq+(ocond-oliq) .EQ. ocond
[3999]	93	! (which is not trivial)
	94	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[4226]	95	!
[4654]	96
	97	! USE de modules contenant des fonctions.
[4651]	98	USE lmdz_cloudth, ONLY : cloudth, cloudth_v3, cloudth_v6, cloudth_mpc
[5007]	99	USE lmdz_lscp_tools, ONLY : calc_qsat_ecmwf, calc_gammasat
	100	USE lmdz_lscp_tools, ONLY : icefrac_lscp, icefrac_lscp_turb
[4664]	101	USE lmdz_lscp_tools, ONLY : fallice_velocity, distance_to_cloud_top
[4380]	102	USE ice_sursat_mod, ONLY : ice_sursat
[4879]	103	USE lmdz_lscp_poprecip, ONLY : poprecip_precld, poprecip_postcld
[3999]	104
[4664]	105	! Use du module lmdz_lscp_ini contenant les constantes
[4666]	106	USE lmdz_lscp_ini, ONLY : prt_level, lunout
[4664]	107	USE lmdz_lscp_ini, ONLY : seuil_neb, niter_lscp, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min
[4910]	108	USE lmdz_lscp_ini, ONLY : ok_radocond_snow, a_tr_sca, cld_expo_con, cld_expo_lsc
[4664]	109	USE lmdz_lscp_ini, ONLY : iflag_cloudth_vert, iflag_rain_incloud_vol, iflag_t_glace, t_glace_min
[4830]	110	USE lmdz_lscp_ini, ONLY : coef_eva, coef_sub,cld_tau_lsc, cld_tau_con, cld_lc_lsc, cld_lc_con
[4664]	111	USE lmdz_lscp_ini, ONLY : iflag_bergeron, iflag_fisrtilp_qsat, iflag_vice, cice_velo, dice_velo
[4910]	112	USE lmdz_lscp_ini, ONLY : iflag_autoconversion, ffallv_con, ffallv_lsc, min_frac_th_cld
[4664]	113	USE lmdz_lscp_ini, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG
[4915]	114	USE lmdz_lscp_ini, ONLY : ok_poprecip
[5007]	115	USE lmdz_lscp_ini, ONLY : iflag_icefrac
[3999]	116	IMPLICIT NONE
	117
	118	!===============================================================================
[4380]	119	! VARIABLES DECLARATION
[3999]	120	!===============================================================================
	121
	122	! INPUT VARIABLES:
	123	!-----------------
	124
[4380]	125	INTEGER, INTENT(IN) :: klon,klev ! number of horizontal grid points and vertical levels
[3999]	126	REAL, INTENT(IN) :: dtime ! time step [s]
[4380]	127	REAL, INTENT(IN) :: missing_val ! missing value for output
[4059]	128
[3999]	129	REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! inter-layer pressure [Pa]
	130	REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! mid-layer pressure [Pa]
[4654]	131	REAL, DIMENSION(klon,klev), INTENT(IN) :: temp ! temperature (K)
[4686]	132	REAL, DIMENSION(klon,klev), INTENT(IN) :: qt ! total specific humidity (in vapor phase in input) [kg/kg]
[5007]	133	REAL, DIMENSION(klon,klev), INTENT(IN) :: qice_save ! ice specific from previous time step [kg/kg]
[3999]	134	INTEGER, INTENT(IN) :: iflag_cld_th ! flag that determines the distribution of convective clouds
	135	INTEGER, INTENT(IN) :: iflag_ice_thermo! flag to activate the ice thermodynamics
[4226]	136	! CR: if iflag_ice_thermo=2, only convection is active
[4072]	137	LOGICAL, INTENT(IN) :: ok_ice_sursat ! flag to determine if ice sursaturation is activated
[4059]	138
[3999]	139	LOGICAL, DIMENSION(klon,klev), INTENT(IN) :: ptconv ! grid points where deep convection scheme is active
	140
	141	!Inputs associated with thermal plumes
[4226]	142
[5007]	143	REAL, DIMENSION(klon,klev), INTENT(IN) :: tv ! virtual potential temperature [K]
	144	REAL, DIMENSION(klon,klev), INTENT(IN) :: qta ! specific humidity within thermals [kg/kg]
	145	REAL, DIMENSION(klon,klev), INTENT(IN) :: fraca ! fraction of thermals within the mesh [-]
	146	REAL, DIMENSION(klon,klev), INTENT(IN) :: pspsk ! exner potential (p/100000)**(R/cp)
	147	REAL, DIMENSION(klon,klev), INTENT(IN) :: tla ! liquid temperature within thermals [K]
	148	REAL, DIMENSION(klon,klev), INTENT(IN) :: tke !--turbulent kinetic energy [m2/s2]
	149	REAL, DIMENSION(klon,klev), INTENT(IN) :: tke_dissip !--TKE dissipation [m2/s3]
[3999]	150
	151	! INPUT/OUTPUT variables
	152	!------------------------
[4562]	153
[5007]	154	REAL, DIMENSION(klon,klev), INTENT(INOUT) :: thl ! liquid potential temperature [K]
	155	REAL, DIMENSION(klon,klev), INTENT(INOUT) :: ratqs ! function of pressure that sets the large-scale
[4059]	156
[4380]	157
[4059]	158	! Input sursaturation en glace
[5007]	159	REAL, DIMENSION(klon,klev), INTENT(INOUT):: rneb_seri ! fraction nuageuse en memoire
[3999]	160
	161	! OUTPUT variables
	162	!-----------------
	163
[5007]	164	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_t ! temperature increment [K]
	165	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_q ! specific humidity increment [kg/kg]
	166	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_ql ! liquid water increment [kg/kg]
	167	REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_qi ! cloud ice mass increment [kg/kg]
	168	REAL, DIMENSION(klon,klev), INTENT(OUT) :: rneb ! cloud fraction [-]
	169	REAL, DIMENSION(klon,klev), INTENT(OUT) :: rneblsvol ! cloud fraction per unit volume [-]
	170	REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfraclr ! precip fraction clear-sky part [-]
	171	REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfracld ! precip fraction cloudy part [-]
	172	REAL, DIMENSION(klon,klev), INTENT(OUT) :: cldfraliq ! liquid fraction of cloud [-]
	173	REAL, DIMENSION(klon,klev), INTENT(OUT) :: sigma2_icefracturb ! Variance of the diagnostic supersaturation distribution (icefrac_turb) [-]
	174	REAL, DIMENSION(klon,klev), INTENT(OUT) :: mean_icefracturb ! Mean of the diagnostic supersaturation distribution (icefrac_turb) [-]
	175	REAL, DIMENSION(klon,klev), INTENT(OUT) :: radocond ! condensed water used in the radiation scheme [kg/kg]
	176	REAL, DIMENSION(klon,klev), INTENT(OUT) :: radicefrac ! ice fraction of condensed water for radiation scheme
	177	REAL, DIMENSION(klon,klev), INTENT(OUT) :: rhcl ! clear-sky relative humidity [-]
	178	REAL, DIMENSION(klon), INTENT(OUT) :: rain ! surface large-scale rainfall [kg/s/m2]
	179	REAL, DIMENSION(klon), INTENT(OUT) :: snow ! surface large-scale snowfall [kg/s/m2]
	180	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsatl ! saturation specific humidity wrt liquid [kg/kg]
	181	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsats ! saturation specific humidity wrt ice [kg/kg]
	182	REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: prfl ! large-scale rainfall flux in the column [kg/s/m2]
	183	REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: psfl ! large-scale snowfall flux in the column [kg/s/m2]
	184	REAL, DIMENSION(klon,klev), INTENT(OUT) :: distcltop ! distance to cloud top [m]
	185	REAL, DIMENSION(klon,klev), INTENT(OUT) :: temp_cltop ! temperature of cloud top [K]
	186	REAL, DIMENSION(klon,klev), INTENT(OUT) :: beta ! conversion rate of condensed water
[4686]	187
[5007]	188	! fraction of aerosol scavenging through impaction and nucleation (for on-line)
[3999]	189
[5007]	190	REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_impa ! scavenging fraction due tu impaction [-]
	191	REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_nucl ! scavenging fraction due tu nucleation [-]
[3999]	192
[4380]	193	! for supersaturation and contrails parameterisation
[3999]	194
[5007]	195	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qclr ! specific total water content in clear sky region [kg/kg]
	196	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcld ! specific total water content in cloudy region [kg/kg]
	197	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qss ! specific total water content in supersat region [kg/kg]
	198	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qvc ! specific vapor content in clouds [kg/kg]
	199	REAL, DIMENSION(klon,klev), INTENT(OUT) :: rnebclr ! mesh fraction of clear sky [-]
	200	REAL, DIMENSION(klon,klev), INTENT(OUT) :: rnebss ! mesh fraction of ISSR [-]
	201	REAL, DIMENSION(klon,klev), INTENT(OUT) :: gamma_ss ! coefficient governing the ice nucleation RHi threshold [-]
	202	REAL, DIMENSION(klon,klev), INTENT(OUT) :: Tcontr ! threshold temperature for contrail formation [K]
	203	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcontr ! threshold humidity for contrail formation [kg/kg]
	204	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcontr2 ! // (2nd expression more consistent with LMDZ expression of q)
	205	REAL, DIMENSION(klon,klev), INTENT(OUT) :: fcontrN ! fraction of grid favourable to non-persistent contrails
	206	REAL, DIMENSION(klon,klev), INTENT(OUT) :: fcontrP ! fraction of grid favourable to persistent contrails
	207	REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sth ! mean saturation deficit in thermals
	208	REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_senv ! mean saturation deficit in environment
	209	REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sigmath ! std of saturation deficit in thermals
	210	REAL, DIMENSION(klon,klev), INTENT(OUT) :: cloudth_sigmaenv ! std of saturation deficit in environment
[3999]	211
[4803]	212	! for POPRECIP
[4686]	213
[4830]	214	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qraindiag !--DIAGNOSTIC specific rain content [kg/kg]
	215	REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsnowdiag !--DIAGNOSTIC specific snow content [kg/kg]
[4819]	216	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqreva !--rain tendendy due to evaporation [kg/kg/s]
	217	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqssub !--snow tendency due to sublimation [kg/kg/s]
	218	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrcol !--rain tendendy due to collection by rain of liquid cloud droplets [kg/kg/s]
	219	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsagg !--snow tendency due to collection of lcoud ice by aggregation [kg/kg/s]
	220	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrauto !--rain tendency due to autoconversion of cloud liquid [kg/kg/s]
	221	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsauto !--snow tendency due to autoconversion of cloud ice [kg/kg/s]
	222	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsrim !--snow tendency due to riming [kg/kg/s]
	223	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsmelt !--snow tendency due to melting [kg/kg/s]
	224	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrmelt !--rain tendency due to melting [kg/kg/s]
	225	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqsfreez !--snow tendency due to freezing [kg/kg/s]
	226	REAL, DIMENSION(klon,klev), INTENT(OUT) :: dqrfreez !--rain tendency due to freezing [kg/kg/s]
[4686]	227
[4803]	228
[3999]	229	! LOCAL VARIABLES:
	230	!----------------
[5007]	231	REAL,DIMENSION(klon) :: qsl, qsi ! saturation threshold at current vertical level
[4654]	232	REAL :: zct, zcl,zexpo
	233	REAL, DIMENSION(klon,klev) :: ctot
	234	REAL, DIMENSION(klon,klev) :: ctot_vol
	235	REAL, DIMENSION(klon) :: zqs, zdqs
	236	REAL :: zdelta, zcor, zcvm5
	237	REAL, DIMENSION(klon) :: zdqsdT_raw
[5007]	238	REAL, DIMENSION(klon) :: gammasat,dgammasatdt ! coefficient to make cold condensation at the correct RH and derivative wrt T
	239	REAL, DIMENSION(klon) :: Tbef,qlbef,DT ! temperature, humidity and temp. variation during lognormal iteration
[4654]	240	REAL :: num,denom
	241	REAL :: cste
[5007]	242	REAL, DIMENSION(klon) :: zpdf_sig,zpdf_k,zpdf_delta ! lognormal parameters
	243	REAL, DIMENSION(klon) :: Zpdf_a,zpdf_b,zpdf_e1,zpdf_e2 ! lognormal intermediate variables
[4654]	244	REAL :: erf
[4910]	245	REAL, DIMENSION(klon) :: zfice_th
[4654]	246	REAL, DIMENSION(klon) :: qcloud, qincloud_mpc
	247	REAL, DIMENSION(klon) :: zrfl, zrfln
	248	REAL :: zqev, zqevt
	249	REAL, DIMENSION(klon) :: zifl, zifln, ziflprev
	250	REAL :: zqev0,zqevi, zqevti
	251	REAL, DIMENSION(klon) :: zoliq, zcond, zq, zqn
	252	REAL, DIMENSION(klon) :: zoliql, zoliqi
	253	REAL, DIMENSION(klon) :: zt
	254	REAL, DIMENSION(klon,klev) :: zrho
	255	REAL, DIMENSION(klon) :: zdz,iwc
	256	REAL :: zchau,zfroi
	257	REAL, DIMENSION(klon) :: zfice,zneb,znebprecip
	258	REAL :: zmelt,zrain,zsnow,zprecip
	259	REAL, DIMENSION(klon) :: dzfice
[5007]	260	REAL, DIMENSION(klon) :: zfice_turb, dzfice_turb
[4654]	261	REAL :: zsolid
	262	REAL, DIMENSION(klon) :: qtot, qzero
	263	REAL, DIMENSION(klon) :: dqsl,dqsi
	264	REAL :: smallestreal
[3999]	265	! Variables for Bergeron process
[4654]	266	REAL :: zcp, coef1, DeltaT, Deltaq, Deltaqprecl
	267	REAL, DIMENSION(klon) :: zqpreci, zqprecl
[3999]	268	! Variables precipitation energy conservation
[4654]	269	REAL, DIMENSION(klon) :: zmqc
	270	REAL :: zalpha_tr
	271	REAL :: zfrac_lessi
	272	REAL, DIMENSION(klon) :: zprec_cond
[4803]	273	REAL :: zmair
[4654]	274	REAL :: zcpair, zcpeau
	275	REAL, DIMENSION(klon) :: zlh_solid
[4803]	276	REAL, DIMENSION(klon) :: ztupnew
[4913]	277	REAL, DIMENSION(klon) :: zqvapclr, zqupnew ! for poprecip evap / subl
[4654]	278	REAL :: zm_solid ! for liquid -> solid conversion
	279	REAL, DIMENSION(klon) :: zrflclr, zrflcld
	280	REAL, DIMENSION(klon) :: d_zrfl_clr_cld, d_zifl_clr_cld
	281	REAL, DIMENSION(klon) :: d_zrfl_cld_clr, d_zifl_cld_clr
	282	REAL, DIMENSION(klon) :: ziflclr, ziflcld
	283	REAL, DIMENSION(klon) :: znebprecipclr, znebprecipcld
	284	REAL, DIMENSION(klon) :: tot_zneb, tot_znebn, d_tot_zneb
	285	REAL, DIMENSION(klon) :: d_znebprecip_clr_cld, d_znebprecip_cld_clr
	286	REAL, DIMENSION(klon,klev) :: velo
	287	REAL :: vr, ffallv
	288	REAL :: qlmpc, qimpc, rnebmpc
	289	REAL, DIMENSION(klon,klev) :: radocondi, radocondl
	290	REAL :: effective_zneb
[5007]	291	REAL, DIMENSION(klon) :: zdistcltop, ztemp_cltop
	292	REAL, DIMENSION(klon) :: zqliq, zqice, zqvapcl ! for icefrac_lscp_turb
[3999]	293
[4654]	294	INTEGER i, k, n, kk, iter
	295	INTEGER, DIMENSION(klon) :: n_i
[4226]	296	INTEGER ncoreczq
[4654]	297	INTEGER, DIMENSION(klon,klev) :: mpc_bl_points
[4803]	298	LOGICAL iftop
[3999]	299
[4654]	300	LOGICAL, DIMENSION(klon) :: lognormale
	301	LOGICAL, DIMENSION(klon) :: keepgoing
[3999]	302
	303	CHARACTER (len = 20) :: modname = 'lscp'
	304	CHARACTER (len = 80) :: abort_message
	305
	306
	307	!===============================================================================
	308	! INITIALISATION
	309	!===============================================================================
	310
[4226]	311	! Few initial checks
[3999]	312
[4654]	313
[5082]	314	IF (iflag_fisrtilp_qsat < 0) THEN
[4226]	315	abort_message = 'lscp cannot be used with iflag_fisrtilp<0'
	316	CALL abort_physic(modname,abort_message,1)
[3999]	317	ENDIF
	318
	319	! Few initialisations
	320
[4226]	321	znebprecip(:)=0.0
[3999]	322	ctot_vol(1:klon,1:klev)=0.0
	323	rneblsvol(1:klon,1:klev)=0.0
[4226]	324	smallestreal=1.e-9
	325	znebprecipclr(:)=0.0
	326	znebprecipcld(:)=0.0
[3999]	327	mpc_bl_points(:,:)=0
	328
	329	IF (prt_level>9) WRITE(lunout,*) 'NUAGES4 A. JAM'
	330
	331	! AA for 'safety' reasons
	332	zalpha_tr = 0.
	333	zfrac_lessi = 0.
[4380]	334	beta(:,:)= 0.
[3999]	335
[4380]	336	! Initialisation of variables:
	337
[3999]	338	prfl(:,:) = 0.0
	339	psfl(:,:) = 0.0
	340	d_t(:,:) = 0.0
	341	d_q(:,:) = 0.0
	342	d_ql(:,:) = 0.0
	343	d_qi(:,:) = 0.0
	344	rneb(:,:) = 0.0
[4530]	345	pfraclr(:,:)=0.0
	346	pfracld(:,:)=0.0
[5007]	347	cldfraliq(:,:)=0.
	348	sigma2_icefracturb(:,:)=0.
	349	mean_icefracturb(:,:)=0.
[4412]	350	radocond(:,:) = 0.0
[3999]	351	radicefrac(:,:) = 0.0
[4226]	352	frac_nucl(:,:) = 1.0
	353	frac_impa(:,:) = 1.0
[3999]	354	rain(:) = 0.0
	355	snow(:) = 0.0
[4114]	356	zoliq(:)=0.0
	357	zfice(:)=0.0
	358	dzfice(:)=0.0
[5007]	359	zfice_turb(:)=0.0
	360	dzfice_turb(:)=0.0
[4114]	361	zqprecl(:)=0.0
	362	zqpreci(:)=0.0
[3999]	363	zrfl(:) = 0.0
	364	zifl(:) = 0.0
[4114]	365	ziflprev(:)=0.0
[3999]	366	zneb(:) = seuil_neb
[4226]	367	zrflclr(:) = 0.0
	368	ziflclr(:) = 0.0
	369	zrflcld(:) = 0.0
	370	ziflcld(:) = 0.0
	371	tot_zneb(:) = 0.0
	372	tot_znebn(:) = 0.0
	373	d_tot_zneb(:) = 0.0
	374	qzero(:) = 0.0
[5007]	375	zdistcltop(:)=0.0
	376	ztemp_cltop(:) = 0.0
[4803]	377	ztupnew(:)=0.0
[4380]	378	!--ice supersaturation
[4226]	379	gamma_ss(:,:) = 1.0
	380	qss(:,:) = 0.0
	381	rnebss(:,:) = 0.0
[4059]	382	Tcontr(:,:) = missing_val
	383	qcontr(:,:) = missing_val
	384	qcontr2(:,:) = missing_val
	385	fcontrN(:,:) = 0.0
	386	fcontrP(:,:) = 0.0
[4654]	387	distcltop(:,:)=0.
	388	temp_cltop(:,:)=0.
[5007]	389
[4803]	390	!-- poprecip
[4830]	391	qraindiag(:,:)= 0.
	392	qsnowdiag(:,:)= 0.
[4819]	393	dqreva(:,:) = 0.
	394	dqrauto(:,:) = 0.
	395	dqrmelt(:,:) = 0.
	396	dqrfreez(:,:) = 0.
	397	dqrcol(:,:) = 0.
	398	dqssub(:,:) = 0.
	399	dqsauto(:,:) = 0.
	400	dqsrim(:,:) = 0.
	401	dqsagg(:,:) = 0.
	402	dqsfreez(:,:) = 0.
	403	dqsmelt(:,:) = 0.
[4913]	404	zqupnew(:) = 0.
	405	zqvapclr(:) = 0.
[3999]	406
[4803]	407
	408
[4392]	409	!c_iso: variable initialisation for iso
	410
	411
[3999]	412	!===============================================================================
	413	! BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM
	414	!===============================================================================
	415
	416	ncoreczq=0
	417
	418	DO k = klev, 1, -1
	419
[5082]	420	IF (k<=klev-1) THEN
[4803]	421	iftop=.false.
	422	ELSE
	423	iftop=.true.
	424	ENDIF
	425
[3999]	426	! Initialisation temperature and specific humidity
[5007]	427	! temp(klon,klev) is not modified by the routine, instead all changes in temperature are made on zt
	428	! at the end of the klon loop, a temperature incremtent d_t due to all processes
	429	! (thermalization, evap/sub incoming precip, cloud formation, precipitation processes) is calculated
	430	! d_t = temperature tendency due to lscp
	431	! The temperature of the overlying layer is updated here because needed for thermalization
[3999]	432	DO i = 1, klon
[4654]	433	zt(i)=temp(i,k)
[4686]	434	zq(i)=qt(i,k)
[4803]	435	IF (.not. iftop) THEN
[4913]	436	ztupnew(i) = temp(i,k+1) + d_t(i,k+1)
	437	zqupnew(i) = qt(i,k+1) + d_q(i,k+1) + d_ql(i,k+1) + d_qi(i,k+1)
	438	!--zqs(i) is the saturation specific humidity in the layer above
	439	zqvapclr(i) = MAX(0., qt(i,k+1) + d_q(i,k+1) - rneb(i,k+1) * zqs(i))
[4803]	440	ENDIF
[4392]	441	!c_iso init of iso
[3999]	442	ENDDO
[4803]	443
	444	!================================================================
	445	! Flag for the new and more microphysical treatment of precipitation from Atelier Nuage (R)
	446	IF (ok_poprecip) THEN
	447
[4879]	448	CALL poprecip_precld(klon, dtime, iftop, paprs(:,k), paprs(:,k+1), pplay(:,k), &
[4803]	449	zt, ztupnew, zq, zmqc, znebprecipclr, znebprecipcld, &
[4913]	450	zqvapclr, zqupnew, &
[4803]	451	zrfl, zrflclr, zrflcld, &
	452	zifl, ziflclr, ziflcld, &
	453	dqreva(:,k),dqssub(:,k) &
	454	)
	455
	456	!================================================================
	457	ELSE
	458
[3999]	459	! --------------------------------------------------------------------
[4803]	460	! P1> Thermalization of precipitation falling from the overlying layer
[3999]	461	! --------------------------------------------------------------------
[4412]	462	! Computes air temperature variation due to enthalpy transported by
[3999]	463	! precipitation. Precipitation is then thermalized with the air in the
	464	! layer.
	465	! The precipitation should remain thermalized throughout the different
[4412]	466	! thermodynamical transformations.
	467	! The corresponding water mass should
[3999]	468	! be added when calculating the layer's enthalpy change with
	469	! temperature
[4412]	470	! See lmdzpedia page todoan
	471	! todoan: check consistency with ice phase
	472	! todoan: understand why several steps
[3999]	473	! ---------------------------------------------------------------------
	474
[4803]	475	IF (iftop) THEN
[3999]	476
	477	DO i = 1, klon
[4803]	478	zmqc(i) = 0.
	479	ENDDO
	480
	481	ELSE
	482
	483	DO i = 1, klon
[3999]	484
[4803]	485	zmair=(paprs(i,k)-paprs(i,k+1))/RG
[3999]	486	! no condensed water so cp=cp(vapor+dry air)
	487	! RVTMP2=rcpv/rcpd-1
	488	zcpair=RCPD(1.0+RVTMP2zq(i))
[4226]	489	zcpeau=RCPD*RVTMP2
	490
[3999]	491	! zmqc: precipitation mass that has to be thermalized with
	492	! layer's air so that precipitation at the ground has the
	493	! same temperature as the lowermost layer
[4803]	494	zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair
[3999]	495	! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
[4803]	496	zt(i) = ( ztupnew(i)zmqc(i)zcpeau + zcpair*zt(i) ) &
[3999]	497	/ (zcpair + zmqc(i)*zcpeau)
[4226]	498
[3999]	499	ENDDO
	500
	501	ENDIF
	502
	503	! --------------------------------------------------------------------
[4803]	504	! P2> Precipitation evaporation/sublimation/melting
[3999]	505	! --------------------------------------------------------------------
[4803]	506	! A part of the precipitation coming from above is evaporated/sublimated/melted.
[3999]	507	! --------------------------------------------------------------------
[4803]	508
[5082]	509	IF (iflag_evap_prec>=1) THEN
[3999]	510
[4072]	511	! Calculation of saturation specific humidity
	512	! depending on temperature:
[4380]	513	CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
[4072]	514	! wrt liquid water
[4380]	515	CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,1,.false.,qsl(:),dqsl(:))
[4072]	516	! wrt ice
[4380]	517	CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,2,.false.,qsi(:),dqsi(:))
[3999]	518
	519	DO i = 1, klon
[4226]	520
[3999]	521	! if precipitation
[5082]	522	IF (zrfl(i)+zifl(i)>0.) THEN
[4226]	523
[4563]	524	! LudoTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec>=4).
[4392]	525	! c_iso: likely important to distinguish cs from neb isotope precipitation
	526
[5082]	527	IF (iflag_evap_prec>=4) THEN
[4226]	528	zrfl(i) = zrflclr(i)
	529	zifl(i) = ziflclr(i)
	530	ENDIF
[3999]	531
[5082]	532	IF (iflag_evap_prec==1) THEN
[3999]	533	znebprecip(i)=zneb(i)
	534	ELSE
	535	znebprecip(i)=MAX(zneb(i),znebprecip(i))
	536	ENDIF
	537
[5082]	538	IF (iflag_evap_prec>4) THEN
[4563]	539	! Max evaporation not to saturate the clear sky precip fraction
	540	! i.e. the fraction where evaporation occurs
	541	zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecipclr(i))
[5082]	542	ELSEIF (iflag_evap_prec == 4) THEN
[4226]	543	! Max evaporation not to saturate the whole mesh
[4563]	544	! Pay attention -> lead to unrealistic and excessive evaporation
[4226]	545	zqev0 = MAX(0.0, zqs(i)-zq(i))
	546	ELSE
	547	! Max evap not to saturate the fraction below the cloud
	548	zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecip(i))
	549	ENDIF
[3999]	550
	551	! Evaporation of liquid precipitation coming from above
	552	! dP/dz=beta(1-q/qsat)sqrt(P)
	553	! formula from Sundquist 1988, Klemp & Wilhemson 1978
[4563]	554	! LTP: evaporation only in the clear sky part (iflag_evap_prec>=4)
[3999]	555
[5082]	556	IF (iflag_evap_prec==3) THEN
[4072]	557	zqevt = znebprecip(i)coef_eva(1.0-zq(i)/qsl(i)) &
[3999]	558	*SQRT(zrfl(i)/max(1.e-4,znebprecip(i))) &
	559	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
[5082]	560	ELSE IF (iflag_evap_prec>=4) THEN
[4072]	561	zqevt = znebprecipclr(i)coef_eva(1.0-zq(i)/qsl(i)) &
[3999]	562	*SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) &
	563	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
	564	ELSE
[4072]	565	zqevt = 1.coef_eva(1.0-zq(i)/qsl(i))*SQRT(zrfl(i)) &
[3999]	566	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
	567	ENDIF
	568
[4226]	569	zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
	570	RGdtime/(paprs(i,k)-paprs(i,k+1))
[3999]	571
	572	! sublimation of the solid precipitation coming from above
[5082]	573	IF (iflag_evap_prec==3) THEN
[4830]	574	zqevti = znebprecip(i)coef_sub(1.0-zq(i)/qsi(i)) &
[3999]	575	*SQRT(zifl(i)/max(1.e-4,znebprecip(i))) &
	576	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
[5082]	577	ELSE IF (iflag_evap_prec>=4) THEN
[4830]	578	zqevti = znebprecipclr(i)coef_sub(1.0-zq(i)/qsi(i)) &
[3999]	579	*SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) &
	580	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
	581	ELSE
[4830]	582	zqevti = 1.coef_sub(1.0-zq(i)/qsi(i))*SQRT(zifl(i)) &
[3999]	583	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
	584	ENDIF
[4226]	585
[3999]	586	zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
[4226]	587	RGdtime/(paprs(i,k)-paprs(i,k+1))
[3999]	588
	589	! A. JAM
	590	! Evaporation limit: we ensure that the layer's fraction below
	591	! the cloud or the whole mesh (depending on iflag_evap_prec)
	592	! does not reach saturation. In this case, we
	593	! redistribute zqev0 conserving the ratio liquid/ice
[4226]	594
[5082]	595	IF (zqevt+zqevti>zqev0) THEN
[3999]	596	zqev=zqev0*zqevt/(zqevt+zqevti)
	597	zqevi=zqev0*zqevti/(zqevt+zqevti)
	598	ELSE
	599	zqev=zqevt
	600	zqevi=zqevti
	601	ENDIF
	602
	603
[4392]	604	! New solid and liquid precipitation fluxes after evap and sublimation
[3999]	605	zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
	606	/RG/dtime)
	607	zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
	608	/RG/dtime)
	609
[4392]	610
[3999]	611	! vapor, temperature, precip fluxes update
[4803]	612	! vapor is updated after evaporation/sublimation (it is increased)
[3999]	613	zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
	614	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
[4803]	615	! zmqc is the total condensed water in the precip flux (it is decreased)
[3999]	616	zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
	617	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
[4803]	618	! air and precip temperature (i.e., gridbox temperature)
	619	! is updated due to latent heat cooling
[3999]	620	zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
	621	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
	622	* RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) &
	623	+ (zifln(i)-zifl(i)) &
	624	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
	625	* RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))
	626
	627	! New values of liquid and solid precipitation
	628	zrfl(i) = zrfln(i)
	629	zifl(i) = zifln(i)
	630
[4392]	631	! c_iso here call_reevap that updates isotopic zrfl, zifl (in inout)
	632	! due to evap + sublim
	633
	634
[5082]	635	IF (iflag_evap_prec>=4) THEN
[4226]	636	zrflclr(i) = zrfl(i)
	637	ziflclr(i) = zifl(i)
[5082]	638	IF(zrflclr(i) + ziflclr(i)<=0) THEN
[4226]	639	znebprecipclr(i) = 0.0
	640	ENDIF
	641	zrfl(i) = zrflclr(i) + zrflcld(i)
	642	zifl(i) = ziflclr(i) + ziflcld(i)
	643	ENDIF
[3999]	644
[4392]	645	! c_iso duplicate for isotopes or loop on isotopes
[3999]	646
[4392]	647	! Melting:
[4412]	648	zmelt = ((zt(i)-RTT)/(ztfondue-RTT)) ! JYG
[3999]	649	! precip fraction that is melted
	650	zmelt = MIN(MAX(zmelt,0.),1.)
[4392]	651
[4412]	652	! update of rainfall and snowfall due to melting
[5082]	653	IF (iflag_evap_prec>=4) THEN
[4226]	654	zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i)
	655	zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i)
	656	zrfl(i)=zrflclr(i)+zrflcld(i)
	657	ELSE
	658	zrfl(i)=zrfl(i)+zmelt*zifl(i)
	659	ENDIF
[4412]	660
	661
[4392]	662	! c_iso: melting of isotopic precipi with zmelt (no fractionation)
[3999]	663
[4803]	664	! Latent heat of melting because of precipitation melting
	665	! NB: the air + precip temperature is simultaneously updated
[3999]	666	zt(i)=zt(i)-zifl(i)zmelt(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
	667	RLMLT/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))
[4869]	668
[5082]	669	IF (iflag_evap_prec>=4) THEN
[4869]	670	ziflclr(i)=ziflclr(i)*(1.-zmelt)
	671	ziflcld(i)=ziflcld(i)*(1.-zmelt)
	672	zifl(i)=ziflclr(i)+ziflcld(i)
	673	ELSE
	674	zifl(i)=zifl(i)*(1.-zmelt)
	675	ENDIF
[3999]	676
	677	ELSE
	678	! if no precip, we reinitialize the cloud fraction used for the precip to 0
	679	znebprecip(i)=0.
	680
	681	ENDIF ! (zrfl(i)+zifl(i).GT.0.)
	682
[4803]	683	ENDDO ! loop on klon
[4226]	684
[3999]	685	ENDIF ! (iflag_evap_prec>=1)
	686
[4803]	687	ENDIF ! (ok_poprecip)
	688
[3999]	689	! --------------------------------------------------------------------
	690	! End precip evaporation
	691	! --------------------------------------------------------------------
	692
	693	! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter
	694	!-------------------------------------------------------
[4226]	695
[4072]	696	qtot(:)=zq(:)+zmqc(:)
[4380]	697	CALL calc_qsat_ecmwf(klon,zt(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
[4072]	698	DO i = 1, klon
[3999]	699	zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
[4059]	700	zdqsdT_raw(i) = zdqs(i)RCPD(1.0+RVTMP2zq(i)) / (RLVTT(1.-zdelta) + RLSTT*zdelta)
[5082]	701	IF (zq(i) < 1.e-15) THEN
[3999]	702	ncoreczq=ncoreczq+1
	703	zq(i)=1.e-15
	704	ENDIF
[4392]	705	! c_iso: do something similar for isotopes
[3999]	706
[4072]	707	ENDDO
[3999]	708
	709	! --------------------------------------------------------------------
	710	! P2> Cloud formation
	711	!---------------------------------------------------------------------
	712	!
	713	! Unlike fisrtilp, we always assume a 'fractional cloud' approach
	714	! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution
	715	! is prescribed and depends on large scale variables and boundary layer
	716	! properties)
	717	! The decrease in condensed part due tu latent heating is taken into
	718	! account
	719	! -------------------------------------------------------------------
	720
	721	! P2.1> With the PDFs (log-normal, bigaussian)
[4424]	722	! cloud properties calculation with the initial values of t and q
[3999]	723	! ----------------------------------------------------------------
	724
	725	! initialise gammasat and qincloud_mpc
	726	gammasat(:)=1.
	727	qincloud_mpc(:)=0.
	728
[5082]	729	IF (iflag_cld_th>=5) THEN
[4910]	730	! Cloud cover and content in meshes affected by shallow convection,
	731	! are retrieved from a bi-gaussian distribution of the saturation deficit
	732	! following Jam et al. 2013
[3999]	733
[5082]	734	IF (iflag_cloudth_vert<=2) THEN
[4910]	735	! Old version of Arnaud Jam
[3999]	736
[4910]	737	CALL cloudth(klon,klev,k,tv, &
	738	zq,qta,fraca, &
	739	qcloud,ctot,pspsk,paprs,pplay,tla,thl, &
[4654]	740	ratqs,zqs,temp, &
[4651]	741	cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
[3999]	742
[4651]	743
[5082]	744	ELSEIF (iflag_cloudth_vert>=3 .AND. iflag_cloudth_vert<=5) THEN
[4910]	745	! Default version of Arnaud Jam
	746
	747	CALL cloudth_v3(klon,klev,k,tv, &
	748	zq,qta,fraca, &
	749	qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
[4654]	750	ratqs,zqs,temp, &
[4651]	751	cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
[3999]	752
	753
[5082]	754	ELSEIF (iflag_cloudth_vert==6) THEN
[4910]	755	! Jean Jouhaud's version, with specific separation between surface and volume
	756	! cloud fraction Decembre 2018
[3999]	757
[4910]	758	CALL cloudth_v6(klon,klev,k,tv, &
	759	zq,qta,fraca, &
	760	qcloud,ctot,ctot_vol,pspsk,paprs,pplay,tla,thl, &
[4654]	761	ratqs,zqs,temp, &
[4651]	762	cloudth_sth,cloudth_senv,cloudth_sigmath,cloudth_sigmaenv)
[3999]	763
[5082]	764	ELSEIF (iflag_cloudth_vert == 7) THEN
[4910]	765	! Updated version of Arnaud Jam (correction by E. Vignon) + adapted treatment
[5007]	766	! for boundary-layer mixed phase clouds
[4910]	767	CALL cloudth_mpc(klon,klev,k,mpc_bl_points,zt,zq,qta(:,k),fraca(:,k), &
	768	pspsk(:,k),paprs(:,k+1),paprs(:,k),pplay(:,k), tla(:,k), &
	769	ratqs(:,k),qcloud,qincloud_mpc,zfice_th,ctot(:,k),ctot_vol(:,k), &
	770	cloudth_sth(:,k),cloudth_senv(:,k),cloudth_sigmath(:,k),cloudth_sigmaenv(:,k))
	771
[3999]	772	ENDIF
	773
	774
	775	DO i=1,klon
	776	rneb(i,k)=ctot(i,k)
	777	rneblsvol(i,k)=ctot_vol(i,k)
	778	zqn(i)=qcloud(i)
	779	ENDDO
	780
	781	ENDIF
	782
[5082]	783	IF (iflag_cld_th <= 4) THEN
[3999]	784
	785	! lognormal
[5007]	786	lognormale(:) = .TRUE.
[3999]	787
[5082]	788	ELSEIF (iflag_cld_th >= 6) THEN
[3999]	789
	790	! lognormal distribution when no thermals
[5007]	791	lognormale(:) = fraca(:,k) < min_frac_th_cld
[3999]	792
	793	ELSE
[4226]	794	! When iflag_cld_th=5, we always assume
[3999]	795	! bi-gaussian distribution
[5007]	796	lognormale(:) = .FALSE.
[3999]	797
	798	ENDIF
	799
	800	DT(:) = 0.
	801	n_i(:)=0
	802	Tbef(:)=zt(:)
	803	qlbef(:)=0.
	804
	805	! Treatment of non-boundary layer clouds (lognormale)
	806	! condensation with qsat(T) variation (adaptation)
[4424]	807	! Iterative resolution to converge towards qsat
	808	! with update of temperature, ice fraction and qsat at
	809	! each iteration
[4226]	810
[4424]	811	! todoan -> sensitivity to iflag_fisrtilp_qsat
[3999]	812	DO iter=1,iflag_fisrtilp_qsat+1
[4226]	813
[3999]	814	DO i=1,klon
	815
[4424]	816	! keepgoing = .true. while convergence is not satisfied
[5082]	817	keepgoing(i)=ABS(DT(i))>DDT0
[4226]	818
[5082]	819	IF ((keepgoing(i) .OR. (n_i(i) == 0)) .AND. lognormale(i)) THEN
[4226]	820
[3999]	821	! if not convergence:
	822
	823	! P2.2.1> cloud fraction and condensed water mass calculation
	824	! Calculated variables:
	825	! rneb : cloud fraction
	826	! zqn : total water within the cloud
	827	! zcond: mean condensed water within the mesh
	828	! rhcl: clear-sky relative humidity
	829	!---------------------------------------------------------------
	830
[4424]	831	! new temperature that only serves in the iteration process:
[3999]	832	Tbef(i)=Tbef(i)+DT(i)
	833
	834	! Rneb, qzn and zcond for lognormal PDFs
[4072]	835	qtot(i)=zq(i)+zmqc(i)
[4226]	836
[4072]	837	ENDIF
[3999]	838
[4072]	839	ENDDO
	840
[4380]	841	! Calculation of saturation specific humidity and ice fraction
[4226]	842	CALL calc_qsat_ecmwf(klon,Tbef(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:))
	843	CALL calc_gammasat(klon,Tbef(:),qtot(:),pplay(:,k),gammasat(:),dgammasatdt(:))
[4072]	844	! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs
	845	zdqs(:) = gammasat(:)zdqs(:)+zqs(:)dgammasatdt(:)
[4562]	846	! cloud phase determination
[5082]	847	IF (iflag_t_glace>=4) THEN
[4562]	848	! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
[5007]	849	CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
[4562]	850	ENDIF
[4072]	851
[5007]	852	CALL icefrac_lscp(klon, zt(:), iflag_ice_thermo, zdistcltop(:),ztemp_cltop(:),zfice(:),dzfice(:))
	853
[4424]	854	DO i=1,klon !todoan : check if loop in i is needed
[4072]	855
[5082]	856	IF ((keepgoing(i) .OR. (n_i(i) == 0)) .AND. lognormale(i)) THEN
[4072]	857
[3999]	858	zpdf_sig(i)=ratqs(i,k)*zq(i)
	859	zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
	860	zpdf_delta(i)=log(zq(i)/(gammasat(i)*zqs(i)))
	861	zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
	862	zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
	863	zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
	864	zpdf_e1(i)=sign(min(ABS(zpdf_e1(i)),5.),zpdf_e1(i))
	865	zpdf_e1(i)=1.-erf(zpdf_e1(i))
	866	zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
	867	zpdf_e2(i)=sign(min(ABS(zpdf_e2(i)),5.),zpdf_e2(i))
	868	zpdf_e2(i)=1.-erf(zpdf_e2(i))
[4226]	869
[5082]	870	IF ((.NOT.ok_ice_sursat).OR.(Tbef(i)>t_glace_min)) THEN
[4059]	871
[5082]	872	IF (zpdf_e1(i)<1.e-10) THEN
[4059]	873	rneb(i,k)=0.
	874	zqn(i)=gammasat(i)*zqs(i)
	875	ELSE
	876	rneb(i,k)=0.5*zpdf_e1(i)
	877	zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
	878	ENDIF
	879
[4072]	880	rnebss(i,k)=0.0 !--added by OB (needed because of the convergence loop on time)
[4059]	881	fcontrN(i,k)=0.0 !--idem
	882	fcontrP(i,k)=0.0 !--idem
	883	qss(i,k)=0.0 !--idem
[4226]	884
[4424]	885	ELSE ! in case of ice supersaturation by Audran
[4226]	886
[4059]	887	!------------------------------------
[4072]	888	! ICE SUPERSATURATION
[4059]	889	!------------------------------------
[3999]	890
[4654]	891	CALL ice_sursat(pplay(i,k), paprs(i,k)-paprs(i,k+1), dtime, i, k, temp(i,k), zq(i), &
[4226]	892	gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k), &
	893	rneb(i,k), zqn(i), rnebss(i,k), qss(i,k), &
[4059]	894	Tcontr(i,k), qcontr(i,k), qcontr2(i,k), fcontrN(i,k), fcontrP(i,k) )
	895
	896	ENDIF ! ((flag_ice_sursat.eq.0).or.(Tbef(i).gt.t_glace_min))
	897
[3999]	898	! If vertical heterogeneity, change fraction by volume as well
[5082]	899	IF (iflag_cloudth_vert>=3) THEN
[3999]	900	ctot_vol(i,k)=rneb(i,k)
	901	rneblsvol(i,k)=ctot_vol(i,k)
	902	ENDIF
	903
	904
[4072]	905	! P2.2.2> Approximative calculation of temperature variation DT
	906	! due to condensation.
	907	! Calculated variables:
	908	! dT : temperature change due to condensation
	909	!---------------------------------------------------------------
[3999]	910
	911
[5082]	912	IF (zfice(i)<1) THEN
[3999]	913	cste=RLVTT
	914	ELSE
	915	cste=RLSTT
	916	ENDIF
[5007]	917
	918	! LEA_R : check formule
[3999]	919	qlbef(i)=max(0.,zqn(i)-zqs(i))
	920	num = -Tbef(i)+zt(i)+rneb(i,k)((1-zfice(i))RLVTT &
[4226]	921	+zfice(i)RLSTT)/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))*qlbef(i)
[3999]	922	denom = 1.+rneb(i,k)((1-zfice(i))RLVTT+zfice(i)RLSTT)/cstezdqs(i) &
	923	-(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2(zq(i)+zmqc(i)))rneb(i,k) &
[4226]	924	qlbef(i)dzfice(i)
[4424]	925	! here we update a provisory temperature variable that only serves in the iteration
	926	! process
[3999]	927	DT(i)=num/denom
	928	n_i(i)=n_i(i)+1
	929
[4424]	930	ENDIF ! end keepgoing
[3999]	931
	932	ENDDO ! end loop on i
	933
	934	ENDDO ! iter=1,iflag_fisrtilp_qsat+1
	935
	936	! P2.3> Final quantities calculation
	937	! Calculated variables:
	938	! rneb : cloud fraction
	939	! zcond: mean condensed water in the mesh
	940	! zqn : mean water vapor in the mesh
[4562]	941	! zfice: ice fraction in clouds
[3999]	942	! zt : temperature
	943	! rhcl : clear-sky relative humidity
	944	! ----------------------------------------------------------------
	945
	946
[4562]	947	! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top
[5082]	948	IF (iflag_t_glace>=4) THEN
[5007]	949	CALL distance_to_cloud_top(klon,klev,k,temp,pplay,paprs,rneb,zdistcltop,ztemp_cltop)
	950	distcltop(:,k)=zdistcltop(:)
	951	temp_cltop(:,k)=ztemp_cltop(:)
	952	ENDIF
[3999]	953
[5007]	954	! Partition function depending on temperature
	955	CALL icefrac_lscp(klon, Tbef, iflag_ice_thermo, zdistcltop, ztemp_cltop, zfice, dzfice)
[4562]	956
[5007]	957	! Partition function depending on tke for non shallow-convective clouds
[5082]	958	IF (iflag_icefrac >= 1) THEN
[5007]	959
	960	CALL icefrac_lscp_turb(klon, dtime, Tbef, pplay(:,k), paprs(:,k), paprs(:,k+1), qice_save(:,k), ziflcld, zqn, &
	961	rneb(:,k), tke(:,k), tke_dissip(:,k), zqliq, zqvapcl, zqice, zfice_turb, dzfice_turb, cldfraliq(:,k),sigma2_icefracturb(:,k), mean_icefracturb(:,k))
	962
	963	ENDIF
	964
[4562]	965	! Water vapor update, Phase determination and subsequent latent heat exchange
[3999]	966	DO i=1, klon
[5007]	967	! Overwrite phase partitioning in boundary layer mixed phase clouds when the
	968	! iflag_cloudth_vert=7 and specific param is activated
[5082]	969	IF (mpc_bl_points(i,k) > 0) THEN
[3999]	970	zcond(i) = MAX(0.0,qincloud_mpc(i))*rneb(i,k)
	971	! following line is very strange and probably wrong
	972	rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i)
	973	! water vapor update and partition function if thermals
	974	zq(i) = zq(i) - zcond(i)
[4910]	975	zfice(i)=zfice_th(i)
[3999]	976	ELSE
	977	! Checks on rneb, rhcl and zqn
[5082]	978	IF (rneb(i,k) <= 0.0) THEN
[3999]	979	zqn(i) = 0.0
	980	rneb(i,k) = 0.0
	981	zcond(i) = 0.0
	982	rhcl(i,k)=zq(i)/zqs(i)
[5082]	983	ELSE IF (rneb(i,k) >= 1.0) THEN
[3999]	984	zqn(i) = zq(i)
	985	rneb(i,k) = 1.0
[4818]	986	zcond(i) = MAX(0.0,zqn(i)-zqs(i))
[3999]	987	rhcl(i,k)=1.0
	988	ELSE
[4818]	989	zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
[3999]	990	! following line is very strange and probably wrong:
	991	rhcl(i,k)=(zqs(i)+zq(i))/2./zqs(i)
[5007]	992	! Overwrite partitioning for non shallow-convective clouds if iflag_icefrac>1 (icefrac turb param)
[5082]	993	IF (iflag_icefrac >= 1) THEN
[5007]	994	IF (lognormale(i)) THEN
	995	zcond(i) = zqliq(i) + zqice(i)
	996	zfice(i)=zfice_turb(i)
	997	rhcl(i,k) = zqvapcl(i) * rneb(i,k) + (zq(i) - zqn(i)) * (1.-rneb(i,k))
	998	ENDIF
	999	ENDIF
[3999]	1000	ENDIF
	1001
[4226]	1002	! water vapor update
	1003	zq(i) = zq(i) - zcond(i)
[3999]	1004
	1005	ENDIF
	1006
[4392]	1007	! c_iso : routine that computes in-cloud supersaturation
	1008	! c_iso condensation of isotopes (zcond, zsursat, zfice, zq in input)
	1009
[3999]	1010	! temperature update due to phase change
	1011	zt(i) = zt(i) + (1.-zfice(i))*zcond(i) &
	1012	& * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) &
	1013	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i)))
	1014	ENDDO
	1015
	1016	! If vertical heterogeneity, change volume fraction
[5082]	1017	IF (iflag_cloudth_vert >= 3) THEN
[3999]	1018	ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.)
	1019	rneblsvol(1:klon,k)=ctot_vol(1:klon,k)
	1020	ENDIF
	1021
	1022	! ----------------------------------------------------------------
[4114]	1023	! P3> Precipitation formation
[3999]	1024	! ----------------------------------------------------------------
[4226]	1025
[4803]	1026	!================================================================
	1027	IF (ok_poprecip) THEN
	1028
	1029	DO i = 1, klon
	1030	zoliql(i) = zcond(i) * ( 1. - zfice(i) )
	1031	zoliqi(i) = zcond(i) * zfice(i)
	1032	ENDDO
	1033
	1034	CALL poprecip_postcld(klon, dtime, paprs(:,k), paprs(:,k+1), pplay(:,k), &
	1035	ctot_vol(:,k), ptconv(:,k), &
	1036	zt, zq, zoliql, zoliqi, zfice, &
	1037	rneb(:,k), znebprecipclr, znebprecipcld, &
	1038	zrfl, zrflclr, zrflcld, &
	1039	zifl, ziflclr, ziflcld, &
[4830]	1040	qraindiag(:,k), qsnowdiag(:,k), dqrauto(:,k), &
[4819]	1041	dqrcol(:,k), dqrmelt(:,k), dqrfreez(:,k), &
	1042	dqsauto(:,k), dqsagg(:,k), dqsrim(:,k), &
	1043	dqsmelt(:,k), dqsfreez(:,k) &
[4803]	1044	)
	1045
	1046	DO i = 1, klon
	1047	zoliq(i) = zoliql(i) + zoliqi(i)
[5082]	1048	IF ( zoliq(i) > 0. ) THEN
[4803]	1049	zfice(i) = zoliqi(i)/zoliq(i)
	1050	ELSE
	1051	zfice(i) = 0.0
	1052	ENDIF
[4819]	1053
	1054	! calculation of specific content of condensates seen by radiative scheme
	1055	IF (ok_radocond_snow) THEN
	1056	radocond(i,k) = zoliq(i)
	1057	radocondl(i,k)= radocond(i,k)*(1.-zfice(i))
[4830]	1058	radocondi(i,k)= radocond(i,k)*zfice(i)+qsnowdiag(i,k)
[4819]	1059	ELSE
	1060	radocond(i,k) = zoliq(i)
	1061	radocondl(i,k)= radocond(i,k)*(1.-zfice(i))
	1062	radocondi(i,k)= radocond(i,k)*zfice(i)
	1063	ENDIF
[4803]	1064	ENDDO
	1065
	1066	!================================================================
	1067	ELSE
	1068
[3999]	1069	! LTP:
[5082]	1070	IF (iflag_evap_prec >= 4) THEN
[3999]	1071
	1072	!Partitionning between precipitation coming from clouds and that coming from CS
	1073
	1074	!0) Calculate tot_zneb, total cloud fraction above the cloud
	1075	!assuming a maximum-random overlap (voir Jakob and Klein, 2000)
	1076
	1077	DO i=1, klon
[4424]	1078	tot_znebn(i) = 1. - (1.-tot_zneb(i))*(1 - max(rneb(i,k),zneb(i))) &
	1079	/(1.-min(zneb(i),1.-smallestreal))
[3999]	1080	d_tot_zneb(i) = tot_znebn(i) - tot_zneb(i)
	1081	tot_zneb(i) = tot_znebn(i)
	1082
	1083
	1084	!1) Cloudy to clear air
	1085	d_znebprecip_cld_clr(i) = znebprecipcld(i) - min(rneb(i,k),znebprecipcld(i))
[5082]	1086	IF (znebprecipcld(i) > 0.) THEN
[3999]	1087	d_zrfl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*zrflcld(i)
	1088	d_zifl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*ziflcld(i)
	1089	ELSE
	1090	d_zrfl_cld_clr(i) = 0.
	1091	d_zifl_cld_clr(i) = 0.
	1092	ENDIF
	1093
	1094	!2) Clear to cloudy air
[4226]	1095	d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i,k) - d_tot_zneb(i) - zneb(i)))
[5082]	1096	IF (znebprecipclr(i) > 0) THEN
[3999]	1097	d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i)
	1098	d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*ziflclr(i)
	1099	ELSE
	1100	d_zrfl_clr_cld(i) = 0.
	1101	d_zifl_clr_cld(i) = 0.
	1102	ENDIF
	1103
	1104	!Update variables
[4226]	1105	znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i)
[3999]	1106	znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i)
	1107	zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i)
	1108	ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i)
	1109	zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i)
	1110	ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i)
	1111
[4392]	1112	! c_iso do the same thing for isotopes precip
[3999]	1113	ENDDO
	1114	ENDIF
	1115
[4559]	1116
	1117	! Autoconversion
	1118	! -------------------------------------------------------------------------------
	1119
	1120
[4412]	1121	! Initialisation of zoliq and radocond variables
[3999]	1122
	1123	DO i = 1, klon
	1124	zoliq(i) = zcond(i)
[4072]	1125	zoliqi(i)= zoliq(i)*zfice(i)
	1126	zoliql(i)= zoliq(i)*(1.-zfice(i))
[4392]	1127	! c_iso : initialisation of zoliq* also for isotopes
[4114]	1128	zrho(i,k) = pplay(i,k) / zt(i) / RD
	1129	zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG)
[4072]	1130	iwc(i) = 0.
	1131	zneb(i) = MAX(rneb(i,k),seuil_neb)
[4559]	1132	radocond(i,k) = zoliq(i)/REAL(niter_lscp+1)
	1133	radocondi(i,k)=zoliq(i)*zfice(i)/REAL(niter_lscp+1)
	1134	radocondl(i,k)=zoliq(i)*(1.-zfice(i))/REAL(niter_lscp+1)
[3999]	1135	ENDDO
	1136
[4559]	1137
	1138	DO n = 1, niter_lscp
[3999]	1139
	1140	DO i=1,klon
[5082]	1141	IF (rneb(i,k)>0.0) THEN
[4114]	1142	iwc(i) = zrho(i,k) * zoliqi(i) / zneb(i) ! in-cloud ice condensate content
[3999]	1143	ENDIF
	1144	ENDDO
	1145
[4226]	1146	CALL fallice_velocity(klon,iwc(:),zt(:),zrho(:,k),pplay(:,k),ptconv(:,k),velo(:,k))
[3999]	1147
	1148	DO i = 1, klon
	1149
[5082]	1150	IF (rneb(i,k)>0.0) THEN
[3999]	1151
[4072]	1152	! Initialization of zrain, zsnow and zprecip:
	1153	zrain=0.
	1154	zsnow=0.
	1155	zprecip=0.
[4392]	1156	! c_iso same init for isotopes. Externalisation?
[3999]	1157
[5082]	1158	IF (zneb(i)==seuil_neb) THEN
[4072]	1159	zprecip = 0.0
	1160	zsnow = 0.0
	1161	zrain= 0.0
[3999]	1162	ELSE
[4420]	1163
[3999]	1164	IF (ptconv(i,k)) THEN ! if convective point
	1165	zcl=cld_lc_con
	1166	zct=1./cld_tau_con
[4559]	1167	zexpo=cld_expo_con
	1168	ffallv=ffallv_con
[3999]	1169	ELSE
	1170	zcl=cld_lc_lsc
	1171	zct=1./cld_tau_lsc
[4559]	1172	zexpo=cld_expo_lsc
	1173	ffallv=ffallv_lsc
[3999]	1174	ENDIF
	1175
[4559]	1176
[3999]	1177	! if vertical heterogeneity is taken into account, we use
	1178	! the "true" volume fraction instead of a modified
	1179	! surface fraction (which is larger and artificially
	1180	! reduces the in-cloud water).
[4072]	1181
[4559]	1182	! Liquid water quantity to remove: zchau (Sundqvist, 1978)
	1183	! dqliq/dt=-qliq/tau(1-exp(-qcin/clw)*2)
	1184	!.........................................................
[5082]	1185	IF ((iflag_cloudth_vert>=3).AND.(iflag_rain_incloud_vol==1)) THEN
[3999]	1186
[4420]	1187	! if vertical heterogeneity is taken into account, we use
	1188	! the "true" volume fraction instead of a modified
	1189	! surface fraction (which is larger and artificially
	1190	! reduces the in-cloud water).
[4559]	1191	effective_zneb=ctot_vol(i,k)
	1192	ELSE
	1193	effective_zneb=zneb(i)
	1194	ENDIF
[4420]	1195
	1196
[5082]	1197	IF (iflag_autoconversion == 2) THEN
[4559]	1198	! two-steps resolution with niter_lscp=1 sufficient
	1199	! we first treat the second term (with exponential) in an explicit way
	1200	! and then treat the first term (-q/tau) in an exact way
	1201	zchau=zoliql(i)(1.-exp(-dtime/REAL(niter_lscp)zct &
	1202	(1.-exp(-(zoliql(i)/effective_zneb/zcl)*zexpo))))
[4420]	1203	ELSE
[4559]	1204	! old explicit resolution with subtimesteps
	1205	zchau = zctdtime/REAL(niter_lscp)zoliql(i) &
	1206	(1.0-EXP(-(zoliql(i)/effective_zneb/zcl)*zexpo))
[4420]	1207	ENDIF
	1208
	1209
[4559]	1210	! Ice water quantity to remove (Zender & Kiehl, 1997)
	1211	! dqice/dt=1/rhod(rhowice*qice)/dz
	1212	!....................................
[5082]	1213	IF (iflag_autoconversion == 2) THEN
[4559]	1214	! exact resolution, niter_lscp=1 is sufficient but works only
	1215	! with iflag_vice=0
[5082]	1216	IF (zoliqi(i) > 0.) THEN
[4559]	1217	zfroi=(zoliqi(i)-((zoliqi(i)**(-dice_velo)) &
	1218	+dice_velodtime/REAL(niter_lscp)cice_velo/zdz(i)ffallv)*(-1./dice_velo))
	1219	ELSE
	1220	zfroi=0.
	1221	ENDIF
	1222	ELSE
	1223	! old explicit resolution with subtimesteps
	1224	zfroi = dtime/REAL(niter_lscp)/zdz(i)zoliqi(i)velo(i,k)
	1225	ENDIF
	1226
[4397]	1227	zrain = MIN(MAX(zchau,0.0),zoliql(i))
	1228	zsnow = MIN(MAX(zfroi,0.0),zoliqi(i))
[4072]	1229	zprecip = MAX(zrain + zsnow,0.0)
[3999]	1230
	1231	ENDIF
[4226]	1232
[5082]	1233	IF (iflag_autoconversion >= 1) THEN
[4559]	1234	! debugged version with phase conservation through the autoconversion process
	1235	zoliql(i) = MAX(zoliql(i)-1.*zrain , 0.0)
	1236	zoliqi(i) = MAX(zoliqi(i)-1.*zsnow , 0.0)
	1237	zoliq(i) = MAX(zoliq(i)-zprecip , 0.0)
	1238	ELSE
	1239	! bugged version with phase resetting
	1240	zoliql(i) = MAX(zoliq(i)(1.-zfice(i))-1.zrain , 0.0)
	1241	zoliqi(i) = MAX(zoliq(i)zfice(i)-1.zsnow , 0.0)
	1242	zoliq(i) = MAX(zoliq(i)-zprecip , 0.0)
	1243	ENDIF
[3999]	1244
[4392]	1245	! c_iso: call isotope_conversion (for readibility) or duplicate
	1246
[4559]	1247	radocond(i,k) = radocond(i,k) + zoliq(i)/REAL(niter_lscp+1)
	1248	radocondl(i,k) = radocondl(i,k) + zoliql(i)/REAL(niter_lscp+1)
	1249	radocondi(i,k) = radocondi(i,k) + zoliqi(i)/REAL(niter_lscp+1)
[4226]	1250
[4072]	1251	ENDIF ! rneb >0
[3999]	1252
	1253	ENDDO ! i = 1,klon
	1254
	1255	ENDDO ! n = 1,niter
	1256
[4072]	1257	! Precipitation flux calculation
[3999]	1258
	1259	DO i = 1, klon
[4380]	1260
[5082]	1261	IF (iflag_evap_prec>=4) THEN
[4380]	1262	ziflprev(i)=ziflcld(i)
	1263	ELSE
	1264	ziflprev(i)=zifl(i)*zneb(i)
	1265	ENDIF
[3999]	1266
[5082]	1267	IF (rneb(i,k) > 0.0) THEN
[3999]	1268
[4072]	1269	! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux:
	1270	! If T<0C, liquid precip are converted into ice, which leads to an increase in
	1271	! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly
	1272	! taken into account through a linearization of the equations and by approximating
	1273	! the liquid precipitation process with a "threshold" process. We assume that
	1274	! condensates are not modified during this operation. Liquid precipitation is
	1275	! removed (in the limit DeltaT<273.15-T). Solid precipitation increases.
	1276	! Water vapor increases as well
	1277	! Note that compared to fisrtilp, we always assume iflag_bergeron=2
	1278
	1279	zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i)
	1280	zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i))
	1281	zcp=RCPD(1.0+RVTMP2(zq(i)+zmqc(i)+zcond(i)))
	1282	coef1 = rneb(i,k)RLSTT/zcpzdqsdT_raw(i)
	1283	! Computation of DT if all the liquid precip freezes
	1284	DeltaT = RLMLTzqprecl(i) / (zcp(1.+coef1))
	1285	! T should not exceed the freezing point
	1286	! that is Delta > RTT-zt(i)
	1287	DeltaT = max( min( RTT-zt(i), DeltaT) , 0. )
	1288	zt(i) = zt(i) + DeltaT
	1289	! water vaporization due to temp. increase
	1290	Deltaq = rneb(i,k)zdqsdT_raw(i)DeltaT
	1291	! we add this vaporized water to the total vapor and we remove it from the precip
	1292	zq(i) = zq(i) + Deltaq
	1293	! The three "max" lines herebelow protect from rounding errors
	1294	zcond(i) = max( zcond(i)- Deltaq, 0. )
	1295	! liquid precipitation converted to ice precip
	1296	Deltaqprecl = -zcp/RLMLT(1.+coef1)DeltaT
	1297	zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. )
	1298	! iced water budget
	1299	zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.)
[4226]	1300
[4392]	1301	! c_iso : mv here condensation of isotopes + redispatchage en precip
	1302
[5082]	1303	IF (iflag_evap_prec>=4) THEN
[4226]	1304	zrflcld(i) = zrflcld(i)+zqprecl(i) &
	1305	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
	1306	ziflcld(i) = ziflcld(i)+ zqpreci(i) &
	1307	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
	1308	znebprecipcld(i) = rneb(i,k)
	1309	zrfl(i) = zrflcld(i) + zrflclr(i)
	1310	zifl(i) = ziflcld(i) + ziflclr(i)
	1311	ELSE
[3999]	1312	zrfl(i) = zrfl(i)+ zqprecl(i) &
	1313	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
	1314	zifl(i) = zifl(i)+ zqpreci(i) &
[4226]	1315	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
[3999]	1316	ENDIF
[4392]	1317	! c_iso : same for isotopes
[3999]	1318
	1319	ENDIF ! rneb>0
	1320
	1321	ENDDO
	1322
[4226]	1323	! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min
[4563]	1324	! if iflag_evap_prec>=4
[5082]	1325	IF (iflag_evap_prec>=4) THEN
[4114]	1326
[4226]	1327	DO i=1,klon
[4114]	1328
[5082]	1329	IF ((zrflclr(i) + ziflclr(i)) > 0. ) THEN
[3999]	1330	znebprecipclr(i) = min(znebprecipclr(i),max(zrflclr(i)/ &
	1331	(MAX(znebprecipclr(i),seuil_neb)rain_int_min), ziflclr(i)/(MAX(znebprecipclr(i),seuil_neb)rain_int_min)))
	1332	ELSE
[4226]	1333	znebprecipclr(i)=0.0
	1334	ENDIF
	1335
[5082]	1336	IF ((zrflcld(i) + ziflcld(i)) > 0.) THEN
[3999]	1337	znebprecipcld(i) = min(znebprecipcld(i), max(zrflcld(i)/ &
	1338	(MAX(znebprecipcld(i),seuil_neb)rain_int_min), ziflcld(i)/(MAX(znebprecipcld(i),seuil_neb)rain_int_min)))
[4226]	1339	ELSE
	1340	znebprecipcld(i)=0.0
	1341	ENDIF
[5007]	1342	!IF ( ((1-zfice(i))*zoliq(i) .GT. 0.) .AND. (zt(i) .LE. 233.15) ) THEN
	1343	!print*,'WARNING LEA OLIQ A <-40°C '
	1344	!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)
	1345	!ENDIF
[4226]	1346	ENDDO
[3999]	1347
[4226]	1348	ENDIF
[3999]	1349
[4910]	1350
[4803]	1351	ENDIF ! ok_poprecip
	1352
[3999]	1353	! End of precipitation formation
	1354	! --------------------------------
	1355
	1356
[4803]	1357	! Calculation of cloud condensates amount seen by radiative scheme
	1358	!-----------------------------------------------------------------
	1359
[4114]	1360	! Calculation of the concentration of condensates seen by the radiation scheme
[4412]	1361	! for liquid phase, we take radocondl
	1362	! for ice phase, we take radocondi if we neglect snowfall, otherwise (ok_radocond_snow=true)
[4803]	1363	! we recalculate radocondi to account for contributions from the precipitation flux
	1364	! TODO: for the moment, we deactivate this option when ok_poprecip=.true.
[3999]	1365
[5082]	1366	IF ((ok_radocond_snow) .AND. (k < klev) .AND. (.NOT. ok_poprecip)) THEN
[4114]	1367	! for the solid phase (crystals + snowflakes)
[4412]	1368	! we recalculate radocondi by summing
[4114]	1369	! the ice content calculated in the mesh
	1370	! + the contribution of the non-evaporated snowfall
	1371	! from the overlying layer
	1372	DO i=1,klon
[5082]	1373	IF (ziflprev(i) /= 0.0) THEN
[4412]	1374	radocondi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)+ziflprev(i)/zrho(i,k+1)/velo(i,k+1)
[4114]	1375	ELSE
[4412]	1376	radocondi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)
[4114]	1377	ENDIF
[4412]	1378	radocond(i,k)=radocondl(i,k)+radocondi(i,k)
[4114]	1379	ENDDO
	1380	ENDIF
	1381
[4412]	1382	! caculate the percentage of ice in "radocond" so cloud+precip seen by the radiation scheme
[4114]	1383	DO i=1,klon
[5082]	1384	IF (radocond(i,k) > 0.) THEN
[4412]	1385	radicefrac(i,k)=MIN(MAX(radocondi(i,k)/radocond(i,k),0.),1.)
[4114]	1386	ENDIF
	1387	ENDDO
	1388
[3999]	1389	! ----------------------------------------------------------------
	1390	! P4> Wet scavenging
	1391	! ----------------------------------------------------------------
	1392
	1393	!Scavenging through nucleation in the layer
	1394
	1395	DO i = 1,klon
	1396
[5082]	1397	IF(zcond(i)>zoliq(i)+1.e-10) THEN
[3999]	1398	beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
	1399	ELSE
	1400	beta(i,k) = 0.
	1401	ENDIF
	1402
[4059]	1403	zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG
[3999]	1404
[5082]	1405	IF (rneb(i,k)>0.0.AND.zprec_cond(i)>0.) THEN
[3999]	1406
[5082]	1407	IF (temp(i,k) >= t_glace_min) THEN
[3999]	1408	zalpha_tr = a_tr_sca(3)
	1409	ELSE
	1410	zalpha_tr = a_tr_sca(4)
	1411	ENDIF
	1412
	1413	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
	1414	frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
[4226]	1415
[3999]	1416	! Nucleation with a factor of -1 instead of -0.5
	1417	zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
	1418
	1419	ENDIF
	1420
[4226]	1421	ENDDO
	1422
[3999]	1423	! Scavenging through impaction in the underlying layer
	1424
	1425	DO kk = k-1, 1, -1
	1426
	1427	DO i = 1, klon
	1428
[5082]	1429	IF (rneb(i,k)>0.0.AND.zprec_cond(i)>0.) THEN
[3999]	1430
[5082]	1431	IF (temp(i,kk) >= t_glace_min) THEN
[3999]	1432	zalpha_tr = a_tr_sca(1)
	1433	ELSE
	1434	zalpha_tr = a_tr_sca(2)
	1435	ENDIF
	1436
	1437	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
	1438	frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
	1439
	1440	ENDIF
	1441
	1442	ENDDO
	1443
	1444	ENDDO
[4226]	1445
[4072]	1446	!--save some variables for ice supersaturation
[4059]	1447	!
	1448	DO i = 1, klon
[4072]	1449	! for memory
[4059]	1450	rneb_seri(i,k) = rneb(i,k)
[3999]	1451
[4072]	1452	! for diagnostics
[4059]	1453	rnebclr(i,k) = 1.0 - rneb(i,k) - rnebss(i,k)
[3999]	1454
[4059]	1455	qvc(i,k) = zqs(i) * rneb(i,k)
[4072]	1456	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
[4059]	1457	qcld(i,k) = qvc(i,k) + zcond(i)
[4226]	1458	ENDDO
[4072]	1459	!q_sat
[4380]	1460	CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs(:))
	1461	CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,2,.false.,qsats(:,k),zdqs(:))
[4059]	1462
[4803]	1463
	1464
	1465	! Outputs:
	1466	!-------------------------------
	1467	! Precipitation fluxes at layer interfaces
	1468	! + precipitation fractions +
	1469	! temperature and water species tendencies
	1470	DO i = 1, klon
	1471	psfl(i,k)=zifl(i)
	1472	prfl(i,k)=zrfl(i)
	1473	pfraclr(i,k)=znebprecipclr(i)
	1474	pfracld(i,k)=znebprecipcld(i)
	1475	d_ql(i,k) = (1-zfice(i))*zoliq(i)
	1476	d_qi(i,k) = zfice(i)*zoliq(i)
	1477	d_q(i,k) = zq(i) - qt(i,k)
	1478	! c_iso: same for isotopes
	1479	d_t(i,k) = zt(i) - temp(i,k)
	1480	ENDDO
	1481
	1482
[4226]	1483	ENDDO
[4059]	1484
[3999]	1485
	1486	! Rain or snow at the surface (depending on the first layer temperature)
	1487	DO i = 1, klon
	1488	snow(i) = zifl(i)
	1489	rain(i) = zrfl(i)
[4392]	1490	! c_iso final output
[3999]	1491	ENDDO
	1492
	1493	IF (ncoreczq>0) THEN
	1494	WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.'
	1495	ENDIF
	1496
[4654]	1497
	1498	RETURN
	1499
[4380]	1500	END SUBROUTINE lscp
[3999]	1501	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	1502
[4664]	1503	END MODULE lmdz_lscp

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format