Context Navigation

lscp_mod.F90 @ 4640

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

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

Download in other formats:

Original Format