Context Navigation

lmdz_lscp.F90 @ 4670

Last change on this file since 4670 was 4666, checked in by fhourdin, 13 months ago
Replayisation lmdz_lscp_old
File size: 55.5 KB

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

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

Download in other formats:

Original Format