Context Navigation

lscp_mod.F90 @ 4553

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

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

Download in other formats:

Original Format