Context Navigation

lscp_mod.F90 @ 4443

Last change on this file since 4443 was 4425, checked in by evignon, 18 months ago
fin des commentaires de lscp a l'issue de l'atelier nuages, noms de variables de contrails (Olivier)
File size: 51.9 KB

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

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

Download in other formats:

Original Format