Context Navigation

lmdz_cloudth.F90 @ 5224

Last change on this file since 5224 was 5160, checked in by abarral, 6 months ago
Put .h into modules
File size: 92.5 KB

Rev	Line
[4651]	1	MODULE lmdz_cloudth
[2686]	2
	3	IMPLICIT NONE
	4
	5	CONTAINS
	6
[5144]	7	SUBROUTINE cloudth(ngrid, klev, ind2, &
	8	ztv, po, zqta, fraca, &
	9	qcloud, ctot, zpspsk, paprs, pplay, ztla, zthl, &
	10	ratqs, zqs, t, &
	11	cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv)
[2686]	12
[5144]	13	USE lmdz_cloudth_ini, ONLY: iflag_cloudth_vert, iflag_ratqs
	14	USE lmdz_yoethf
[5153]	15
[5144]	16	USE lmdz_yomcst
[2686]	17
[5144]	18	IMPLICIT NONE
[5153]	19	INCLUDE "FCTTRE.h"
[4535]	20
[2686]	21
[5144]	22	!===========================================================================
	23	! Auteur : Arnaud Octavio Jam (LMD/CNRS)
	24	! Date : 25 Mai 2010
	25	! Objet : calcule les valeurs de qc et rneb dans les thermiques
	26	!===========================================================================
[2686]	27
[5144]	28	INTEGER itap, ind1, ind2
	29	INTEGER ngrid, klev, klon, l, ig
	30	REAL, DIMENSION(ngrid, klev), INTENT(OUT) :: cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv
[2686]	31
[5144]	32	REAL ztv(ngrid, klev)
	33	REAL po(ngrid)
	34	REAL zqenv(ngrid)
	35	REAL zqta(ngrid, klev)
[2686]	36
[5144]	37	REAL fraca(ngrid, klev + 1)
	38	REAL zpspsk(ngrid, klev)
	39	REAL paprs(ngrid, klev + 1)
	40	REAL pplay(ngrid, klev)
	41	REAL ztla(ngrid, klev)
	42	REAL zthl(ngrid, klev)
[2686]	43
[5144]	44	REAL zqsatth(ngrid, klev)
	45	REAL zqsatenv(ngrid, klev)
[2686]	46
[5144]	47	REAL sigma1(ngrid, klev)
	48	REAL sigma2(ngrid, klev)
	49	REAL qlth(ngrid, klev)
	50	REAL qlenv(ngrid, klev)
	51	REAL qltot(ngrid, klev)
	52	REAL cth(ngrid, klev)
	53	REAL cenv(ngrid, klev)
	54	REAL ctot(ngrid, klev)
	55	REAL rneb(ngrid, klev)
	56	REAL t(ngrid, klev)
	57	REAL qsatmmussig1, qsatmmussig2, sqrt2pi, pi
	58	REAL rdd, cppd, Lv
	59	REAL alth, alenv, ath, aenv
	60	REAL sth, senv, sigma1s, sigma2s, xth, xenv
	61	REAL Tbef, zdelta, qsatbef, zcor
	62	REAL qlbef
	63	REAL ratqs(ngrid, klev) ! determine la largeur de distribution de vapeur
[2686]	64
[5144]	65	REAL zpdf_sig(ngrid), zpdf_k(ngrid), zpdf_delta(ngrid)
	66	REAL zpdf_a(ngrid), zpdf_b(ngrid), zpdf_e1(ngrid), zpdf_e2(ngrid)
	67	REAL zqs(ngrid), qcloud(ngrid)
	68	REAL erf
[2686]	69
	70
	71
	72
[5144]	73	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	74	! Gestion de deux versions de cloudth
	75	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	76
	77	IF (iflag_cloudth_vert>=1) THEN
	78	CALL cloudth_vert(ngrid, klev, ind2, &
	79	ztv, po, zqta, fraca, &
	80	qcloud, ctot, zpspsk, paprs, pplay, ztla, zthl, &
	81	ratqs, zqs, t)
[2686]	82	RETURN
[5144]	83	ENDIF
	84	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[2686]	85
	86
[5144]	87	!-------------------------------------------------------------------------------
	88	! Initialisation des variables r?elles
	89	!-------------------------------------------------------------------------------
	90	sigma1(:, :) = 0.
	91	sigma2(:, :) = 0.
	92	qlth(:, :) = 0.
	93	qlenv(:, :) = 0.
	94	qltot(:, :) = 0.
	95	rneb(:, :) = 0.
	96	qcloud(:) = 0.
	97	cth(:, :) = 0.
	98	cenv(:, :) = 0.
	99	ctot(:, :) = 0.
	100	qsatmmussig1 = 0.
	101	qsatmmussig2 = 0.
	102	rdd = 287.04
	103	cppd = 1005.7
	104	pi = 3.14159
	105	Lv = 2.5e6
	106	sqrt2pi = sqrt(2. * pi)
[2686]	107
	108
	109
[5144]	110	!-------------------------------------------------------------------------------
	111	! Calcul de la fraction du thermique et des ?cart-types des distributions
	112	!-------------------------------------------------------------------------------
[5158]	113	DO ind1 = 1, ngrid
[2686]	114
[5144]	115	IF ((ztv(ind1, 1)>ztv(ind1, 2)).AND.(fraca(ind1, ind2)>1.e-10)) THEN
	116	zqenv(ind1) = (po(ind1) - fraca(ind1, ind2) * zqta(ind1, ind2)) / (1. - fraca(ind1, ind2))
[2686]	117
	118
[5144]	119	! zqenv(ind1)=po(ind1)
	120	Tbef = zthl(ind1, ind2) * zpspsk(ind1, ind2)
	121	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	122	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	123	qsatbef = MIN(0.5, qsatbef)
	124	zcor = 1. / (1. - retv * qsatbef)
	125	qsatbef = qsatbef * zcor
	126	zqsatenv(ind1, ind2) = qsatbef
[2686]	127
[5144]	128	alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2)
	129	aenv = 1. / (1. + (alenv * Lv / cppd))
	130	senv = aenv * (po(ind1) - zqsatenv(ind1, ind2))
[2686]	131
[5144]	132	Tbef = ztla(ind1, ind2) * zpspsk(ind1, ind2)
	133	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	134	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	135	qsatbef = MIN(0.5, qsatbef)
	136	zcor = 1. / (1. - retv * qsatbef)
	137	qsatbef = qsatbef * zcor
	138	zqsatth(ind1, ind2) = qsatbef
[2686]	139
[5144]	140	alth = (0.622 * Lv * zqsatth(ind1, ind2)) / (rdd * ztla(ind1, ind2)**2)
	141	ath = 1. / (1. + (alth * Lv / cppd))
	142	sth = ath * (zqta(ind1, ind2) - zqsatth(ind1, ind2))
[2686]	143
	144
	145
[5144]	146	!------------------------------------------------------------------------------
	147	! Calcul des ?cart-types pour s
	148	!------------------------------------------------------------------------------
[2686]	149
[5144]	150	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
	151	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)0.4+0.002zqta(ind1,ind2)
	152	! if (paprs(ind1,ind2).gt.90000) THEN
	153	! ratqs(ind1,ind2)=0.002
	154	! else
	155	! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
	156	! endif
	157	sigma1s = (1.1*0.5) (fraca(ind1, ind2)*0.6) / (1 - fraca(ind1, ind2)) ((sth - senv)2)0.5 + 0.002 * po(ind1)
	158	sigma2s = 0.11 * ((sth - senv)2)0.5 / (fraca(ind1, ind2) + 0.01)*0.4 + 0.002 zqta(ind1, ind2)
	159	! sigma1s=ratqs(ind1,ind2)*po(ind1)
	160	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)*0.4+0.00003
[2686]	161
[5144]	162	!------------------------------------------------------------------------------
	163	! Calcul de l'eau condens?e et de la couverture nuageuse
	164	!------------------------------------------------------------------------------
	165	sqrt2pi = sqrt(2. * pi)
	166	xth = sth / (sqrt(2.) * sigma2s)
	167	xenv = senv / (sqrt(2.) * sigma1s)
	168	cth(ind1, ind2) = 0.5 * (1. + 1. * erf(xth))
	169	cenv(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv))
	170	ctot(ind1, ind2) = fraca(ind1, ind2) * cth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv(ind1, ind2)
[2686]	171
[5144]	172	qlth(ind1, ind2) = sigma2s * ((exp(-1. * xth*2) / sqrt2pi) + xth sqrt(2.) * cth(ind1, ind2))
	173	qlenv(ind1, ind2) = sigma1s * ((exp(-1. * xenv*2) / sqrt2pi) + xenv sqrt(2.) * cenv(ind1, ind2))
	174	qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2)
[2686]	175
[5144]	176	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	177	IF (ctot(ind1, ind2)<1.e-10) THEN
	178	ctot(ind1, ind2) = 0.
	179	qcloud(ind1) = zqsatenv(ind1, ind2)
[2686]	180
[5144]	181	else
[2686]	182
[5144]	183	ctot(ind1, ind2) = ctot(ind1, ind2)
	184	qcloud(ind1) = qltot(ind1, ind2) / ctot(ind1, ind2) + zqs(ind1)
[2686]	185
[5144]	186	endif
[2686]	187
	188
[5144]	189	! PRINT*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
[2686]	190
	191	else ! gaussienne environnement seule
	192
[5144]	193	zqenv(ind1) = po(ind1)
	194	Tbef = t(ind1, ind2)
	195	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	196	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	197	qsatbef = MIN(0.5, qsatbef)
	198	zcor = 1. / (1. - retv * qsatbef)
	199	qsatbef = qsatbef * zcor
	200	zqsatenv(ind1, ind2) = qsatbef
[2686]	201
	202
[5144]	203	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
	204	zthl(ind1, ind2) = t(ind1, ind2) * (101325 / paprs(ind1, ind2))**(rdd / cppd)
	205	alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2)
	206	aenv = 1. / (1. + (alenv * Lv / cppd))
	207	senv = aenv * (po(ind1) - zqsatenv(ind1, ind2))
[2686]	208
[5144]	209	sigma1s = ratqs(ind1, ind2) * zqenv(ind1)
[2686]	210
[5144]	211	sqrt2pi = sqrt(2. * pi)
	212	xenv = senv / (sqrt(2.) * sigma1s)
	213	ctot(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv))
	214	qltot(ind1, ind2) = sigma1s * ((exp(-1. * xenv*2) / sqrt2pi) + xenv sqrt(2.) * cenv(ind1, ind2))
[2686]	215
[5144]	216	IF (ctot(ind1, ind2)<1.e-3) THEN
	217	ctot(ind1, ind2) = 0.
	218	qcloud(ind1) = zqsatenv(ind1, ind2)
[2686]	219
[5144]	220	else
[2686]	221
[5144]	222	ctot(ind1, ind2) = ctot(ind1, ind2)
	223	qcloud(ind1) = qltot(ind1, ind2) / ctot(ind1, ind2) + zqsatenv(ind1, ind2)
[2686]	224
[5144]	225	endif
[2686]	226
[5144]	227	endif
	228	enddo
[2686]	229
[5144]	230	RETURN
	231	! end
	232	END SUBROUTINE cloudth
[2686]	233
	234
[5144]	235	!===========================================================================
	236	SUBROUTINE cloudth_vert(ngrid, klev, ind2, &
	237	ztv, po, zqta, fraca, &
	238	qcloud, ctot, zpspsk, paprs, pplay, ztla, zthl, &
	239	ratqs, zqs, t)
[5143]	240
[5144]	241	!===========================================================================
	242	! Auteur : Arnaud Octavio Jam (LMD/CNRS)
	243	! Date : 25 Mai 2010
	244	! Objet : calcule les valeurs de qc et rneb dans les thermiques
	245	!===========================================================================
[2686]	246
[5144]	247	USE lmdz_cloudth_ini, ONLY: iflag_cloudth_vert, vert_alpha
	248	USE lmdz_yoethf
[5153]	249
[5144]	250	USE lmdz_yomcst
[2686]	251
[5144]	252	IMPLICIT NONE
[5153]	253	INCLUDE "FCTTRE.h"
[2686]	254
[5144]	255	INTEGER itap, ind1, ind2
	256	INTEGER ngrid, klev, klon, l, ig
[2686]	257
[5144]	258	REAL ztv(ngrid, klev)
	259	REAL po(ngrid)
	260	REAL zqenv(ngrid)
	261	REAL zqta(ngrid, klev)
[2686]	262
[5144]	263	REAL fraca(ngrid, klev + 1)
	264	REAL zpspsk(ngrid, klev)
	265	REAL paprs(ngrid, klev + 1)
	266	REAL pplay(ngrid, klev)
	267	REAL ztla(ngrid, klev)
	268	REAL zthl(ngrid, klev)
[2686]	269
[5144]	270	REAL zqsatth(ngrid, klev)
	271	REAL zqsatenv(ngrid, klev)
[2686]	272
[5144]	273	REAL sigma1(ngrid, klev)
	274	REAL sigma2(ngrid, klev)
	275	REAL qlth(ngrid, klev)
	276	REAL qlenv(ngrid, klev)
	277	REAL qltot(ngrid, klev)
	278	REAL cth(ngrid, klev)
	279	REAL cenv(ngrid, klev)
	280	REAL ctot(ngrid, klev)
	281	REAL rneb(ngrid, klev)
	282	REAL t(ngrid, klev)
	283	REAL qsatmmussig1, qsatmmussig2, sqrt2pi, pi
	284	REAL rdd, cppd, Lv, sqrt2, sqrtpi
	285	REAL alth, alenv, ath, aenv
	286	REAL sth, senv, sigma1s, sigma2s, xth, xenv
	287	REAL xth1, xth2, xenv1, xenv2, deltasth, deltasenv
	288	REAL IntJ, IntI1, IntI2, IntI3, coeffqlenv, coeffqlth
	289	REAL Tbef, zdelta, qsatbef, zcor
	290	REAL qlbef
	291	REAL ratqs(ngrid, klev) ! determine la largeur de distribution de vapeur
	292	! Change the width of the PDF used for vertical subgrid scale heterogeneity
	293	! (J Jouhaud, JL Dufresne, JB Madeleine)
[2686]	294
[5144]	295	REAL zpdf_sig(ngrid), zpdf_k(ngrid), zpdf_delta(ngrid)
	296	REAL zpdf_a(ngrid), zpdf_b(ngrid), zpdf_e1(ngrid), zpdf_e2(ngrid)
	297	REAL zqs(ngrid), qcloud(ngrid)
	298	REAL erf
[2686]	299
[5144]	300	!------------------------------------------------------------------------------
	301	! Initialisation des variables r?elles
	302	!------------------------------------------------------------------------------
	303	sigma1(:, :) = 0.
	304	sigma2(:, :) = 0.
	305	qlth(:, :) = 0.
	306	qlenv(:, :) = 0.
	307	qltot(:, :) = 0.
	308	rneb(:, :) = 0.
	309	qcloud(:) = 0.
	310	cth(:, :) = 0.
	311	cenv(:, :) = 0.
	312	ctot(:, :) = 0.
	313	qsatmmussig1 = 0.
	314	qsatmmussig2 = 0.
	315	rdd = 287.04
	316	cppd = 1005.7
	317	pi = 3.14159
	318	Lv = 2.5e6
	319	sqrt2pi = sqrt(2. * pi)
	320	sqrt2 = sqrt(2.)
	321	sqrtpi = sqrt(pi)
[2686]	322
[5144]	323	!-------------------------------------------------------------------------------
	324	! Calcul de la fraction du thermique et des ?cart-types des distributions
	325	!-------------------------------------------------------------------------------
[5158]	326	DO ind1 = 1, ngrid
[2686]	327
[5144]	328	IF ((ztv(ind1, 1)>ztv(ind1, 2)).AND.(fraca(ind1, ind2)>1.e-10)) THEN
	329	zqenv(ind1) = (po(ind1) - fraca(ind1, ind2) * zqta(ind1, ind2)) / (1. - fraca(ind1, ind2))
[2686]	330
	331
[5144]	332	! zqenv(ind1)=po(ind1)
	333	Tbef = zthl(ind1, ind2) * zpspsk(ind1, ind2)
	334	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	335	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	336	qsatbef = MIN(0.5, qsatbef)
	337	zcor = 1. / (1. - retv * qsatbef)
	338	qsatbef = qsatbef * zcor
	339	zqsatenv(ind1, ind2) = qsatbef
[2686]	340
[5144]	341	alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2)
	342	aenv = 1. / (1. + (alenv * Lv / cppd))
	343	senv = aenv * (po(ind1) - zqsatenv(ind1, ind2))
[2686]	344
[5144]	345	Tbef = ztla(ind1, ind2) * zpspsk(ind1, ind2)
	346	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	347	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	348	qsatbef = MIN(0.5, qsatbef)
	349	zcor = 1. / (1. - retv * qsatbef)
	350	qsatbef = qsatbef * zcor
	351	zqsatth(ind1, ind2) = qsatbef
[2686]	352
[5144]	353	alth = (0.622 * Lv * zqsatth(ind1, ind2)) / (rdd * ztla(ind1, ind2)**2)
	354	ath = 1. / (1. + (alth * Lv / cppd))
	355	sth = ath * (zqta(ind1, ind2) - zqsatth(ind1, ind2))
[2686]	356
	357
	358
[5144]	359	!------------------------------------------------------------------------------
	360	! Calcul des ?cart-types pour s
	361	!------------------------------------------------------------------------------
[2686]	362
[5144]	363	sigma1s = (0.92*0.5) (fraca(ind1, ind2)*0.5) / (1 - fraca(ind1, ind2)) ((sth - senv)2)0.5 + ratqs(ind1, ind2) * po(ind1)
	364	sigma2s = 0.09 * ((sth - senv)2)0.5 / (fraca(ind1, ind2) + 0.02)*0.5 + 0.002 zqta(ind1, ind2)
	365	! if (paprs(ind1,ind2).gt.90000) THEN
	366	! ratqs(ind1,ind2)=0.002
	367	! else
	368	! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
	369	! endif
	370	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
	371	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
	372	! sigma1s=ratqs(ind1,ind2)*po(ind1)
	373	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)*0.4+0.00003
[2686]	374
[5144]	375	!------------------------------------------------------------------------------
	376	! Calcul de l'eau condens?e et de la couverture nuageuse
	377	!------------------------------------------------------------------------------
	378	sqrt2pi = sqrt(2. * pi)
	379	xth = sth / (sqrt(2.) * sigma2s)
	380	xenv = senv / (sqrt(2.) * sigma1s)
	381	cth(ind1, ind2) = 0.5 * (1. + 1. * erf(xth))
	382	cenv(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv))
	383	ctot(ind1, ind2) = fraca(ind1, ind2) * cth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv(ind1, ind2)
[2686]	384
[5144]	385	qlth(ind1, ind2) = sigma2s * ((exp(-1. * xth*2) / sqrt2pi) + xth sqrt(2.) * cth(ind1, ind2))
	386	qlenv(ind1, ind2) = sigma1s * ((exp(-1. * xenv*2) / sqrt2pi) + xenv sqrt(2.) * cenv(ind1, ind2))
	387	qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2)
[2686]	388
[5144]	389	IF (iflag_cloudth_vert == 1) THEN
	390	!-------------------------------------------------------------------------------
	391	! Version 2: Modification selon J.-Louis. On condense ?? partir de qsat-ratqs
	392	!-------------------------------------------------------------------------------
	393	! deltasenv=aenvratqs(ind1,ind2)po(ind1)
	394	! deltasth=athratqs(ind1,ind2)zqta(ind1,ind2)
	395	deltasenv = aenv * ratqs(ind1, ind2) * zqsatenv(ind1, ind2)
	396	deltasth = ath * ratqs(ind1, ind2) * zqsatth(ind1, ind2)
	397	! deltasenv=aenv0.01po(ind1)
	398	! deltasth=ath0.01zqta(ind1,ind2)
	399	xenv1 = (senv - deltasenv) / (sqrt(2.) * sigma1s)
	400	xenv2 = (senv + deltasenv) / (sqrt(2.) * sigma1s)
	401	xth1 = (sth - deltasth) / (sqrt(2.) * sigma2s)
	402	xth2 = (sth + deltasth) / (sqrt(2.) * sigma2s)
	403	coeffqlenv = (sigma1s)*2 / (2 sqrtpi * deltasenv)
	404	coeffqlth = (sigma2s)*2 / (2 sqrtpi * deltasth)
[2686]	405
[5144]	406	cth(ind1, ind2) = 0.5 * (1. + 1. * erf(xth2))
	407	cenv(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv2))
	408	ctot(ind1, ind2) = fraca(ind1, ind2) * cth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv(ind1, ind2)
[2686]	409
[5144]	410	IntJ = sigma1s * (exp(-1. * xenv1*2) / sqrt2pi) + 0.5 senv * (1 + erf(xenv1))
	411	IntI1 = coeffqlenv * 0.5 * (0.5 * sqrtpi * (erf(xenv2) - erf(xenv1)) + xenv1 * exp(-1. * xenv1*2) - xenv2 exp(-1. * xenv2**2))
	412	IntI2 = coeffqlenv * xenv2 * (exp(-1. * xenv2*2) - exp(-1. xenv1**2))
	413	IntI3 = coeffqlenv * 0.5 * sqrtpi * xenv2*2 (erf(xenv2) - erf(xenv1))
[2686]	414
[5144]	415	qlenv(ind1, ind2) = IntJ + IntI1 + IntI2 + IntI3
	416	! qlenv(ind1,ind2)=IntJ
	417	! PRINT*, qlenv(ind1,ind2),'VERIF EAU'
[2686]	418
[5144]	419	IntJ = sigma2s * (exp(-1. * xth1*2) / sqrt2pi) + 0.5 sth * (1 + erf(xth1))
	420	! IntI1=coeffqlth((0.5xth1-xth2)exp(-1.xth1*2)+0.5xth2exp(-1.xth2**2))
	421	! IntI2=coeffqlth0.5sqrtpi(0.5+xth22)(erf(xth2)-erf(xth1))
	422	IntI1 = coeffqlth * 0.5 * (0.5 * sqrtpi * (erf(xth2) - erf(xth1)) + xth1 * exp(-1. * xth1*2) - xth2 exp(-1. * xth2**2))
	423	IntI2 = coeffqlth * xth2 * (exp(-1. * xth2*2) - exp(-1. xth1**2))
	424	IntI3 = coeffqlth * 0.5 * sqrtpi * xth2*2 (erf(xth2) - erf(xth1))
	425	qlth(ind1, ind2) = IntJ + IntI1 + IntI2 + IntI3
	426	! qlth(ind1,ind2)=IntJ
	427	! PRINT*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
	428	qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2)
[2686]	429
[5144]	430	ELSE IF (iflag_cloudth_vert == 2) THEN
[2686]	431
[5144]	432	!-------------------------------------------------------------------------------
	433	! Version 3: Modification Jean Jouhaud. On condense a partir de -delta s
	434	!-------------------------------------------------------------------------------
	435	! deltasenv=aenvratqs(ind1,ind2)po(ind1)
	436	! deltasth=athratqs(ind1,ind2)zqta(ind1,ind2)
	437	! deltasenv=aenvratqs(ind1,ind2)zqsatenv(ind1,ind2)
	438	! deltasth=athratqs(ind1,ind2)zqsatth(ind1,ind2)
	439	deltasenv = aenv * vert_alpha * sigma1s
	440	deltasth = ath * vert_alpha * sigma2s
[2686]	441
[5144]	442	xenv1 = -(senv + deltasenv) / (sqrt(2.) * sigma1s)
	443	xenv2 = -(senv - deltasenv) / (sqrt(2.) * sigma1s)
	444	xth1 = -(sth + deltasth) / (sqrt(2.) * sigma2s)
	445	xth2 = -(sth - deltasth) / (sqrt(2.) * sigma2s)
	446	! coeffqlenv=(sigma1s)*2/(2sqrtpi*deltasenv)
	447	! coeffqlth=(sigma2s)*2/(2sqrtpi*deltasth)
[2686]	448
[5144]	449	cth(ind1, ind2) = 0.5 * (1. - 1. * erf(xth1))
	450	cenv(ind1, ind2) = 0.5 * (1. - 1. * erf(xenv1))
	451	ctot(ind1, ind2) = fraca(ind1, ind2) * cth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv(ind1, ind2)
[2686]	452
[5144]	453	IntJ = 0.5 * senv * (1 - erf(xenv2)) + (sigma1s / sqrt2pi) * exp(-1. * xenv2**2)
	454	IntI1 = (((senv + deltasenv)2 + (sigma1s)2) / (8 * deltasenv)) * (erf(xenv2) - erf(xenv1))
	455	IntI2 = (sigma1s*2 / (4 deltasenv * sqrtpi)) * (xenv1 * exp(-1. * xenv1*2) - xenv2 exp(-1. * xenv2**2))
	456	IntI3 = ((sqrt2 * sigma1s * (senv + deltasenv)) / (4 * sqrtpi * deltasenv)) * (exp(-1. * xenv1*2) - exp(-1. xenv2**2))
[2686]	457
[5144]	458	! IntI1=0.5(0.5sqrtpi(erf(xenv2)-erf(xenv1))+xenv1exp(-1.xenv12)-xenv2exp(-1.xenv2*2))
	459	! IntI2=xenv2(exp(-1.xenv2*2)-exp(-1.xenv1**2))
	460	! IntI3=0.5sqrtpixenv2*2(erf(xenv2)-erf(xenv1))
[2686]	461
[5144]	462	qlenv(ind1, ind2) = IntJ + IntI1 + IntI2 + IntI3
	463	! qlenv(ind1,ind2)=IntJ
	464	! PRINT*, qlenv(ind1,ind2),'VERIF EAU'
[2686]	465
[5144]	466	IntJ = 0.5 * sth * (1 - erf(xth2)) + (sigma2s / sqrt2pi) * exp(-1. * xth2**2)
	467	IntI1 = (((sth + deltasth)2 + (sigma2s)2) / (8 * deltasth)) * (erf(xth2) - erf(xth1))
	468	IntI2 = (sigma2s*2 / (4 deltasth * sqrtpi)) * (xth1 * exp(-1. * xth1*2) - xth2 exp(-1. * xth2**2))
	469	IntI3 = ((sqrt2 * sigma2s * (sth + deltasth)) / (4 * sqrtpi * deltasth)) * (exp(-1. * xth1*2) - exp(-1. xth2**2))
[2686]	470
[5144]	471	qlth(ind1, ind2) = IntJ + IntI1 + IntI2 + IntI3
	472	! qlth(ind1,ind2)=IntJ
	473	! PRINT*, IntJ,IntI1,IntI2,IntI3,qlth(ind1,ind2),'VERIF EAU2'
	474	qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2)
	475
	476	ENDIF ! of if (iflag_cloudth_vert==1 or 2)
	477
	478	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	479
	480	IF (cenv(ind1, ind2)<1.e-10.OR.cth(ind1, ind2)<1.e-10) THEN
	481	ctot(ind1, ind2) = 0.
	482	qcloud(ind1) = zqsatenv(ind1, ind2)
	483
	484	else
	485
	486	ctot(ind1, ind2) = ctot(ind1, ind2)
	487	qcloud(ind1) = qltot(ind1, ind2) / ctot(ind1, ind2) + zqs(ind1)
	488	! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
	489	! & +(1.-1.fraca(ind1,ind2))qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
	490
	491	endif
	492
	493
	494
	495	! PRINT*,sth,sigma2s,qlth(ind1,ind2),ctot(ind1,ind2),qltot(ind1,ind2),'verif'
	496
[2686]	497	else ! gaussienne environnement seule
	498
[5144]	499	zqenv(ind1) = po(ind1)
	500	Tbef = t(ind1, ind2)
	501	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	502	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	503	qsatbef = MIN(0.5, qsatbef)
	504	zcor = 1. / (1. - retv * qsatbef)
	505	qsatbef = qsatbef * zcor
	506	zqsatenv(ind1, ind2) = qsatbef
[2686]	507
	508
[5144]	509	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
	510	zthl(ind1, ind2) = t(ind1, ind2) * (101325 / paprs(ind1, ind2))**(rdd / cppd)
	511	alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2)
	512	aenv = 1. / (1. + (alenv * Lv / cppd))
	513	senv = aenv * (po(ind1) - zqsatenv(ind1, ind2))
[2686]	514
[5144]	515	sigma1s = ratqs(ind1, ind2) * zqenv(ind1)
[2686]	516
[5144]	517	sqrt2pi = sqrt(2. * pi)
	518	xenv = senv / (sqrt(2.) * sigma1s)
	519	ctot(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv))
	520	qltot(ind1, ind2) = sigma1s * ((exp(-1. * xenv*2) / sqrt2pi) + xenv sqrt(2.) * cenv(ind1, ind2))
[2686]	521
[5144]	522	IF (ctot(ind1, ind2)<1.e-3) THEN
	523	ctot(ind1, ind2) = 0.
	524	qcloud(ind1) = zqsatenv(ind1, ind2)
[3493]	525
[5144]	526	else
[3493]	527
[5144]	528	ctot(ind1, ind2) = ctot(ind1, ind2)
	529	qcloud(ind1) = qltot(ind1, ind2) / ctot(ind1, ind2) + zqsatenv(ind1, ind2)
[3493]	530
[5144]	531	endif
[2686]	532
[5144]	533	endif
	534	enddo
[2686]	535
[5144]	536	RETURN
	537	! end
	538	END SUBROUTINE cloudth_vert
[2686]	539
	540
[5144]	541	SUBROUTINE cloudth_v3(ngrid, klev, ind2, &
	542	ztv, po, zqta, fraca, &
	543	qcloud, ctot, ctot_vol, zpspsk, paprs, pplay, ztla, zthl, &
	544	ratqs, zqs, t, &
	545	cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv)
[2686]	546
[5144]	547	USE lmdz_cloudth_ini, ONLY: iflag_cloudth_vert
	548	USE lmdz_yoethf
[5153]	549
[5144]	550	USE lmdz_yomcst
[2686]	551
[5144]	552	IMPLICIT NONE
[5153]	553	INCLUDE "FCTTRE.h"
[2686]	554
[4674]	555
[5144]	556	!===========================================================================
	557	! Author : Arnaud Octavio Jam (LMD/CNRS)
	558	! Date : 25 Mai 2010
	559	! Objet : calcule les valeurs de qc et rneb dans les thermiques
	560	!===========================================================================
[4674]	561
[5144]	562	INTEGER, INTENT(IN) :: ind2
	563	INTEGER, INTENT(IN) :: ngrid, klev
[2686]	564
[5144]	565	REAL, DIMENSION(ngrid, klev), INTENT(IN) :: ztv
	566	REAL, DIMENSION(ngrid), INTENT(IN) :: po
	567	REAL, DIMENSION(ngrid, klev), INTENT(IN) :: zqta
	568	REAL, DIMENSION(ngrid, klev + 1), INTENT(IN) :: fraca
	569	REAL, DIMENSION(ngrid), INTENT(OUT) :: qcloud
	570	REAL, DIMENSION(ngrid, klev), INTENT(OUT) :: ctot
	571	REAL, DIMENSION(ngrid, klev), INTENT(OUT) :: ctot_vol
	572	REAL, DIMENSION(ngrid, klev), INTENT(IN) :: zpspsk
	573	REAL, DIMENSION(ngrid, klev + 1), INTENT(IN) :: paprs
	574	REAL, DIMENSION(ngrid, klev), INTENT(IN) :: pplay
	575	REAL, DIMENSION(ngrid, klev), INTENT(IN) :: ztla
	576	REAL, DIMENSION(ngrid, klev), INTENT(INOUT) :: zthl
	577	REAL, DIMENSION(ngrid, klev), INTENT(IN) :: ratqs
	578	REAL, DIMENSION(ngrid), INTENT(IN) :: zqs
	579	REAL, DIMENSION(ngrid, klev), INTENT(IN) :: t
	580	REAL, DIMENSION(ngrid, klev), INTENT(OUT) :: cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv
[2686]	581
[5144]	582	REAL zqenv(ngrid)
	583	REAL zqsatth(ngrid, klev)
	584	REAL zqsatenv(ngrid, klev)
[4674]	585
[5144]	586	REAL sigma1(ngrid, klev)
	587	REAL sigma2(ngrid, klev)
	588	REAL qlth(ngrid, klev)
	589	REAL qlenv(ngrid, klev)
	590	REAL qltot(ngrid, klev)
	591	REAL cth(ngrid, klev)
	592	REAL cenv(ngrid, klev)
	593	REAL cth_vol(ngrid, klev)
	594	REAL cenv_vol(ngrid, klev)
	595	REAL rneb(ngrid, klev)
	596	REAL qsatmmussig1, qsatmmussig2, sqrt2pi, sqrt2, sqrtpi, pi
	597	REAL rdd, cppd, Lv
	598	REAL alth, alenv, ath, aenv
	599	REAL sth, senv, sigma1s, sigma2s, xth, xenv, exp_xenv1, exp_xenv2, exp_xth1, exp_xth2
	600	REAL inverse_rho, beta, a_Brooks, b_Brooks, A_Maj_Brooks, Dx_Brooks, f_Brooks
	601	REAL Tbef, zdelta, qsatbef, zcor
	602	REAL qlbef
	603	REAL zpdf_sig(ngrid), zpdf_k(ngrid), zpdf_delta(ngrid)
	604	REAL zpdf_a(ngrid), zpdf_b(ngrid), zpdf_e1(ngrid), zpdf_e2(ngrid)
	605	REAL erf
	606
	607	INTEGER :: ind1, l, ig
	608
	609	IF (iflag_cloudth_vert>=1) THEN
	610	CALL cloudth_vert_v3(ngrid, klev, ind2, &
	611	ztv, po, zqta, fraca, &
	612	qcloud, ctot, ctot_vol, zpspsk, paprs, pplay, ztla, zthl, &
	613	ratqs, zqs, t, &
	614	cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv)
[2686]	615	RETURN
[5144]	616	ENDIF
	617	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
[2686]	618
	619
[5144]	620	!-------------------------------------------------------------------------------
	621	! Initialisation des variables r?elles
	622	!-------------------------------------------------------------------------------
	623	sigma1(:, :) = 0.
	624	sigma2(:, :) = 0.
	625	qlth(:, :) = 0.
	626	qlenv(:, :) = 0.
	627	qltot(:, :) = 0.
	628	rneb(:, :) = 0.
	629	qcloud(:) = 0.
	630	cth(:, :) = 0.
	631	cenv(:, :) = 0.
	632	ctot(:, :) = 0.
	633	cth_vol(:, :) = 0.
	634	cenv_vol(:, :) = 0.
	635	ctot_vol(:, :) = 0.
	636	qsatmmussig1 = 0.
	637	qsatmmussig2 = 0.
	638	rdd = 287.04
	639	cppd = 1005.7
	640	pi = 3.14159
	641	Lv = 2.5e6
	642	sqrt2pi = sqrt(2. * pi)
	643	sqrt2 = sqrt(2.)
	644	sqrtpi = sqrt(pi)
[2686]	645
	646
[5144]	647	!-------------------------------------------------------------------------------
	648	! Cloud fraction in the thermals and standard deviation of the PDFs
	649	!-------------------------------------------------------------------------------
[5158]	650	DO ind1 = 1, ngrid
[2686]	651
[5144]	652	IF ((ztv(ind1, 1)>ztv(ind1, 2)).AND.(fraca(ind1, ind2)>1.e-10)) THEN
	653	zqenv(ind1) = (po(ind1) - fraca(ind1, ind2) * zqta(ind1, ind2)) / (1. - fraca(ind1, ind2))
[2686]	654
[5144]	655	Tbef = zthl(ind1, ind2) * zpspsk(ind1, ind2)
	656	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	657	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	658	qsatbef = MIN(0.5, qsatbef)
	659	zcor = 1. / (1. - retv * qsatbef)
	660	qsatbef = qsatbef * zcor
	661	zqsatenv(ind1, ind2) = qsatbef
[2686]	662
[5144]	663	alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2) !qsl, p84
	664	aenv = 1. / (1. + (alenv * Lv / cppd)) !al, p84
	665	senv = aenv * (po(ind1) - zqsatenv(ind1, ind2)) !s, p84
[2686]	666
[5144]	667	!po = qt de l'environnement ET des thermique
	668	!zqenv = qt environnement
	669	!zqsatenv = qsat environnement
	670	!zthl = Tl environnement
[2686]	671
[5144]	672	Tbef = ztla(ind1, ind2) * zpspsk(ind1, ind2)
	673	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	674	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	675	qsatbef = MIN(0.5, qsatbef)
	676	zcor = 1. / (1. - retv * qsatbef)
	677	qsatbef = qsatbef * zcor
	678	zqsatth(ind1, ind2) = qsatbef
[2686]	679
[5144]	680	alth = (0.622 * Lv * zqsatth(ind1, ind2)) / (rdd * ztla(ind1, ind2)**2) !qsl, p84
	681	ath = 1. / (1. + (alth * Lv / cppd)) !al, p84
	682	sth = ath * (zqta(ind1, ind2) - zqsatth(ind1, ind2)) !s, p84
[2686]	683
[5144]	684	!zqta = qt thermals
	685	!zqsatth = qsat thermals
	686	!ztla = Tl thermals
[2686]	687
[5144]	688	!------------------------------------------------------------------------------
	689	! s standard deviations
	690	!------------------------------------------------------------------------------
[2686]	691
[5144]	692	! tests
	693	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
	694	! sigma1s=(0.92(fraca(ind1,ind2)0.5)/(1-fraca(ind1,ind2))(((sth-senv)2)0.5))+ratqs(ind1,ind2)*po(ind1)
	695	! sigma2s=(0.09(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)0.5))+0.002zqta(ind1,ind2)
	696	! final option
	697	sigma1s = (1.1*0.5) (fraca(ind1, ind2)*0.6) / (1 - fraca(ind1, ind2)) ((sth - senv)2)0.5 + ratqs(ind1, ind2) * po(ind1)
	698	sigma2s = 0.11 * ((sth - senv)2)0.5 / (fraca(ind1, ind2) + 0.01)*0.4 + 0.002 zqta(ind1, ind2)
[2686]	699
[5144]	700	!------------------------------------------------------------------------------
	701	! Condensed water and cloud cover
	702	!------------------------------------------------------------------------------
	703	xth = sth / (sqrt2 * sigma2s)
	704	xenv = senv / (sqrt2 * sigma1s)
	705	cth(ind1, ind2) = 0.5 * (1. + 1. * erf(xth)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
	706	cenv(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv)) !4.18 p 111, l.7 p115 & 4.20 p 119 thesis Arnaud Jam
	707	ctot(ind1, ind2) = fraca(ind1, ind2) * cth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv(ind1, ind2)
	708	ctot_vol(ind1, ind2) = ctot(ind1, ind2)
[2686]	709
[5144]	710	qlth(ind1, ind2) = sigma2s * ((exp(-1. * xth*2) / sqrt2pi) + xth sqrt2 * cth(ind1, ind2))
	711	qlenv(ind1, ind2) = sigma1s * ((exp(-1. * xenv*2) / sqrt2pi) + xenv sqrt2 * cenv(ind1, ind2))
	712	qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2)
[2686]	713
[5144]	714	IF (ctot(ind1, ind2)<1.e-10) THEN
	715	ctot(ind1, ind2) = 0.
	716	qcloud(ind1) = zqsatenv(ind1, ind2)
	717	else
	718	qcloud(ind1) = qltot(ind1, ind2) / ctot(ind1, ind2) + zqs(ind1)
	719	endif
	720
[2686]	721	else ! Environnement only, follow the if l.110
	722
[5144]	723	zqenv(ind1) = po(ind1)
	724	Tbef = t(ind1, ind2)
	725	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	726	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	727	qsatbef = MIN(0.5, qsatbef)
	728	zcor = 1. / (1. - retv * qsatbef)
	729	qsatbef = qsatbef * zcor
	730	zqsatenv(ind1, ind2) = qsatbef
[2686]	731
[5144]	732	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
	733	zthl(ind1, ind2) = t(ind1, ind2) * (101325 / paprs(ind1, ind2))**(rdd / cppd)
	734	alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2)
	735	aenv = 1. / (1. + (alenv * Lv / cppd))
	736	senv = aenv * (po(ind1) - zqsatenv(ind1, ind2))
[2686]	737
[5144]	738	sigma1s = ratqs(ind1, ind2) * zqenv(ind1)
[2686]	739
[5144]	740	xenv = senv / (sqrt2 * sigma1s)
	741	ctot(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv))
	742	ctot_vol(ind1, ind2) = ctot(ind1, ind2)
	743	qltot(ind1, ind2) = sigma1s * ((exp(-1. * xenv*2) / sqrt2pi) + xenv sqrt2 * cenv(ind1, ind2))
[2686]	744
[5144]	745	IF (ctot(ind1, ind2)<1.e-3) THEN
	746	ctot(ind1, ind2) = 0.
	747	qcloud(ind1) = zqsatenv(ind1, ind2)
	748	else
	749	qcloud(ind1) = qltot(ind1, ind2) / ctot(ind1, ind2) + zqsatenv(ind1, ind2)
	750	endif
	751
[2686]	752	endif ! From the separation (thermal/envrionnement) et (environnement) only, l.110 et l.183
[5144]	753	enddo ! from the loop on ngrid l.108
	754	RETURN
	755	! end
	756	END SUBROUTINE cloudth_v3
[2686]	757
	758
[5144]	759	!===========================================================================
	760	SUBROUTINE cloudth_vert_v3(ngrid, klev, ind2, &
	761	ztv, po, zqta, fraca, &
	762	qcloud, ctot, ctot_vol, zpspsk, paprs, pplay, ztla, zthl, &
	763	ratqs, zqs, t, &
	764	cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv)
[2686]	765
[5144]	766	!===========================================================================
	767	! Auteur : Arnaud Octavio Jam (LMD/CNRS)
	768	! Date : 25 Mai 2010
	769	! Objet : calcule les valeurs de qc et rneb dans les thermiques
	770	!===========================================================================
[2686]	771
[5144]	772	USE lmdz_cloudth_ini, ONLY: iflag_cloudth_vert, iflag_ratqs
	773	USE lmdz_cloudth_ini, ONLY: vert_alpha, vert_alpha_th, sigma1s_factor, sigma1s_power, sigma2s_factor, sigma2s_power, cloudth_ratqsmin, iflag_cloudth_vert_noratqs
	774	USE lmdz_yoethf
[5153]	775
[5144]	776	USE lmdz_yomcst
[2686]	777
[5144]	778	IMPLICIT NONE
[5153]	779	INCLUDE "FCTTRE.h"
[2686]	780
[5144]	781	INTEGER itap, ind1, ind2
	782	INTEGER ngrid, klev, klon, l, ig
	783	REAL, DIMENSION(ngrid, klev), INTENT(OUT) :: cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv
[2686]	784
[5144]	785	REAL ztv(ngrid, klev)
	786	REAL po(ngrid)
	787	REAL zqenv(ngrid)
	788	REAL zqta(ngrid, klev)
[4651]	789
[5144]	790	REAL fraca(ngrid, klev + 1)
	791	REAL zpspsk(ngrid, klev)
	792	REAL paprs(ngrid, klev + 1)
	793	REAL pplay(ngrid, klev)
	794	REAL ztla(ngrid, klev)
	795	REAL zthl(ngrid, klev)
[2686]	796
[5144]	797	REAL zqsatth(ngrid, klev)
	798	REAL zqsatenv(ngrid, klev)
[2957]	799
[5144]	800	REAL sigma1(ngrid, klev)
	801	REAL sigma2(ngrid, klev)
	802	REAL qlth(ngrid, klev)
	803	REAL qlenv(ngrid, klev)
	804	REAL qltot(ngrid, klev)
	805	REAL cth(ngrid, klev)
	806	REAL cenv(ngrid, klev)
	807	REAL ctot(ngrid, klev)
	808	REAL cth_vol(ngrid, klev)
	809	REAL cenv_vol(ngrid, klev)
	810	REAL ctot_vol(ngrid, klev)
	811	REAL rneb(ngrid, klev)
	812	REAL t(ngrid, klev)
	813	REAL qsatmmussig1, qsatmmussig2, sqrtpi, sqrt2, sqrt2pi, pi
	814	REAL rdd, cppd, Lv
	815	REAL alth, alenv, ath, aenv
	816	REAL sth, senv, sigma1s, sigma2s, sigma1s_fraca, sigma1s_ratqs
	817	REAL inverse_rho, beta, a_Brooks, b_Brooks, A_Maj_Brooks, Dx_Brooks, f_Brooks
	818	REAL xth, xenv, exp_xenv1, exp_xenv2, exp_xth1, exp_xth2
	819	REAL xth1, xth2, xenv1, xenv2, deltasth, deltasenv
	820	REAL IntJ, IntI1, IntI2, IntI3, IntJ_CF, IntI1_CF, IntI3_CF, coeffqlenv, coeffqlth
	821	REAL Tbef, zdelta, qsatbef, zcor
	822	REAL qlbef
	823	REAL ratqs(ngrid, klev) ! determine la largeur de distribution de vapeur
	824	! Change the width of the PDF used for vertical subgrid scale heterogeneity
	825	! (J Jouhaud, JL Dufresne, JB Madeleine)
[2686]	826
[5144]	827	REAL zpdf_sig(ngrid), zpdf_k(ngrid), zpdf_delta(ngrid)
	828	REAL zpdf_a(ngrid), zpdf_b(ngrid), zpdf_e1(ngrid), zpdf_e2(ngrid)
	829	REAL zqs(ngrid), qcloud(ngrid)
	830	REAL erf
[3493]	831
[5144]	832	REAL rhodz(ngrid, klev)
	833	REAL zrho(ngrid, klev)
	834	REAL dz(ngrid, klev)
[3493]	835
[5144]	836	DO ind1 = 1, ngrid
	837	!Layer calculation
	838	rhodz(ind1, ind2) = (paprs(ind1, ind2) - paprs(ind1, ind2 + 1)) / rg !kg/m2
	839	zrho(ind1, ind2) = pplay(ind1, ind2) / t(ind1, ind2) / rd !kg/m3
	840	dz(ind1, ind2) = rhodz(ind1, ind2) / zrho(ind1, ind2) !m : epaisseur de la couche en metre
	841	END DO
[3999]	842
[5144]	843	!------------------------------------------------------------------------------
	844	! Initialize
	845	!------------------------------------------------------------------------------
[2686]	846
[5144]	847	sigma1(:, :) = 0.
	848	sigma2(:, :) = 0.
	849	qlth(:, :) = 0.
	850	qlenv(:, :) = 0.
	851	qltot(:, :) = 0.
	852	rneb(:, :) = 0.
	853	qcloud(:) = 0.
	854	cth(:, :) = 0.
	855	cenv(:, :) = 0.
	856	ctot(:, :) = 0.
	857	cth_vol(:, :) = 0.
	858	cenv_vol(:, :) = 0.
	859	ctot_vol(:, :) = 0.
	860	qsatmmussig1 = 0.
	861	qsatmmussig2 = 0.
	862	rdd = 287.04
	863	cppd = 1005.7
	864	pi = 3.14159
	865	Lv = 2.5e6
	866	sqrt2pi = sqrt(2. * pi)
	867	sqrt2 = sqrt(2.)
	868	sqrtpi = sqrt(pi)
[2957]	869
[2686]	870
	871
[5144]	872	!-------------------------------------------------------------------------------
	873	! Calcul de la fraction du thermique et des ecart-types des distributions
	874	!-------------------------------------------------------------------------------
[5158]	875	DO ind1 = 1, ngrid
[2686]	876
[5144]	877	IF ((ztv(ind1, 1)>ztv(ind1, 2)).AND.(fraca(ind1, ind2)>1.e-10)) then !Thermal and environnement
[2686]	878
[5144]	879	zqenv(ind1) = (po(ind1) - fraca(ind1, ind2) * zqta(ind1, ind2)) / (1. - fraca(ind1, ind2)) !qt = aqtth + (1-a)qtenv
[2686]	880
[5144]	881	Tbef = zthl(ind1, ind2) * zpspsk(ind1, ind2)
	882	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	883	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	884	qsatbef = MIN(0.5, qsatbef)
	885	zcor = 1. / (1. - retv * qsatbef)
	886	qsatbef = qsatbef * zcor
	887	zqsatenv(ind1, ind2) = qsatbef
[2686]	888
[5144]	889	alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2) !qsl, p84
	890	aenv = 1. / (1. + (alenv * Lv / cppd)) !al, p84
	891	senv = aenv * (po(ind1) - zqsatenv(ind1, ind2)) !s, p84
[2686]	892
[5144]	893	!zqenv = qt environnement
	894	!zqsatenv = qsat environnement
	895	!zthl = Tl environnement
[2686]	896
[5144]	897	Tbef = ztla(ind1, ind2) * zpspsk(ind1, ind2)
	898	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	899	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	900	qsatbef = MIN(0.5, qsatbef)
	901	zcor = 1. / (1. - retv * qsatbef)
	902	qsatbef = qsatbef * zcor
	903	zqsatth(ind1, ind2) = qsatbef
[2686]	904
[5144]	905	alth = (0.622 * Lv * zqsatth(ind1, ind2)) / (rdd * ztla(ind1, ind2)**2) !qsl, p84
	906	ath = 1. / (1. + (alth * Lv / cppd)) !al, p84
	907	sth = ath * (zqta(ind1, ind2) - zqsatth(ind1, ind2)) !s, p84
[2686]	908
	909
[5144]	910	!zqta = qt thermals
	911	!zqsatth = qsat thermals
	912	!ztla = Tl thermals
	913	!------------------------------------------------------------------------------
	914	! s standard deviation
	915	!------------------------------------------------------------------------------
[2686]	916
[5144]	917	sigma1s_fraca = (sigma1s_factor*0.5) (fraca(ind1, ind2)**sigma1s_power) / &
	918	(1 - fraca(ind1, ind2)) * ((sth - senv)2)0.5
	919	! sigma1s_fraca = (1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5
	920	IF (cloudth_ratqsmin>0.) THEN
	921	sigma1s_ratqs = cloudth_ratqsmin * po(ind1)
	922	ELSE
	923	sigma1s_ratqs = ratqs(ind1, ind2) * po(ind1)
	924	ENDIF
	925	sigma1s = sigma1s_fraca + sigma1s_ratqs
	926	IF (iflag_ratqs==11) THEN
	927	sigma1s = ratqs(ind1, ind2) * po(ind1) * aenv
	928	ENDIF
	929	sigma2s = (sigma2s_factor * (((sth - senv)2)0.5) / ((fraca(ind1, ind2) + 0.02)*sigma2s_power)) + 0.002 zqta(ind1, ind2)
	930	! tests
	931	! sigma1s=(0.92*0.5)(fraca(ind1,ind2)*0.5)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+ratqs(ind1,ind2)*po(ind1)
	932	! sigma1s=(0.92(fraca(ind1,ind2)0.5)/(1-fraca(ind1,ind2))(((sth-senv)2)0.5))+0.002*zqenv(ind1)
	933	! sigma2s=0.09((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)0.5+0.002zqta(ind1,ind2)
	934	! sigma2s=(0.09(((sth-senv)2)0.5)/((fraca(ind1,ind2)+0.02)0.5))+ratqs(ind1,ind2)zqta(ind1,ind2)
	935	! if (paprs(ind1,ind2).gt.90000) THEN
	936	! ratqs(ind1,ind2)=0.002
	937	! else
	938	! ratqs(ind1,ind2)=0.002+0.0*(90000-paprs(ind1,ind2))/20000
	939	! endif
	940	! sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
	941	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
	942	! sigma1s=ratqs(ind1,ind2)*po(ind1)
	943	! sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.02)*0.4+0.00003
[2686]	944
[5144]	945	IF (iflag_cloudth_vert == 1) THEN
	946	!-------------------------------------------------------------------------------
	947	! Version 2: Modification from Arnaud Jam according to JL Dufrense. Condensate from qsat-ratqs
	948	!-------------------------------------------------------------------------------
[2686]	949
[5144]	950	deltasenv = aenv * ratqs(ind1, ind2) * zqsatenv(ind1, ind2)
	951	deltasth = ath * ratqs(ind1, ind2) * zqsatth(ind1, ind2)
[2686]	952
[5144]	953	xenv1 = (senv - deltasenv) / (sqrt(2.) * sigma1s)
	954	xenv2 = (senv + deltasenv) / (sqrt(2.) * sigma1s)
	955	xth1 = (sth - deltasth) / (sqrt(2.) * sigma2s)
	956	xth2 = (sth + deltasth) / (sqrt(2.) * sigma2s)
	957	coeffqlenv = (sigma1s)*2 / (2 sqrtpi * deltasenv)
	958	coeffqlth = (sigma2s)*2 / (2 sqrtpi * deltasth)
[2686]	959
[5144]	960	cth(ind1, ind2) = 0.5 * (1. + 1. * erf(xth2))
	961	cenv(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv2))
	962	ctot(ind1, ind2) = fraca(ind1, ind2) * cth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv(ind1, ind2)
[2686]	963
[5144]	964	! Environment
	965	IntJ = sigma1s * (exp(-1. * xenv1*2) / sqrt2pi) + 0.5 senv * (1 + erf(xenv1))
	966	IntI1 = coeffqlenv * 0.5 * (0.5 * sqrtpi * (erf(xenv2) - erf(xenv1)) + xenv1 * exp(-1. * xenv1*2) - xenv2 exp(-1. * xenv2**2))
	967	IntI2 = coeffqlenv * xenv2 * (exp(-1. * xenv2*2) - exp(-1. xenv1**2))
	968	IntI3 = coeffqlenv * 0.5 * sqrtpi * xenv2*2 (erf(xenv2) - erf(xenv1))
[2686]	969
[5144]	970	qlenv(ind1, ind2) = IntJ + IntI1 + IntI2 + IntI3
[2945]	971
[5144]	972	! Thermal
	973	IntJ = sigma2s * (exp(-1. * xth1*2) / sqrt2pi) + 0.5 sth * (1 + erf(xth1))
	974	IntI1 = coeffqlth * 0.5 * (0.5 * sqrtpi * (erf(xth2) - erf(xth1)) + xth1 * exp(-1. * xth1*2) - xth2 exp(-1. * xth2**2))
	975	IntI2 = coeffqlth * xth2 * (exp(-1. * xth2*2) - exp(-1. xth1**2))
	976	IntI3 = coeffqlth * 0.5 * sqrtpi * xth2*2 (erf(xth2) - erf(xth1))
	977	qlth(ind1, ind2) = IntJ + IntI1 + IntI2 + IntI3
	978	qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2)
[2686]	979
[5144]	980	ELSE IF (iflag_cloudth_vert >= 3) THEN
	981	IF (iflag_cloudth_vert < 5) THEN
	982	!-------------------------------------------------------------------------------
	983	! Version 3: Changes by J. Jouhaud; condensation for q > -delta s
	984	!-------------------------------------------------------------------------------
	985	! deltasenv=aenvratqs(ind1,ind2)po(ind1)
	986	! deltasth=athratqs(ind1,ind2)zqta(ind1,ind2)
	987	! deltasenv=aenvratqs(ind1,ind2)zqsatenv(ind1,ind2)
	988	! deltasth=athratqs(ind1,ind2)zqsatth(ind1,ind2)
	989	IF (iflag_cloudth_vert == 3) THEN
	990	deltasenv = aenv * vert_alpha * sigma1s
	991	deltasth = ath * vert_alpha_th * sigma2s
	992	ELSE IF (iflag_cloudth_vert == 4) THEN
	993	IF (iflag_cloudth_vert_noratqs == 1) THEN
	994	deltasenv = vert_alpha * max(sigma1s_fraca, 1e-10)
	995	deltasth = vert_alpha_th * sigma2s
	996	ELSE
	997	deltasenv = vert_alpha * sigma1s
	998	deltasth = vert_alpha_th * sigma2s
	999	ENDIF
	1000	ENDIF
[2686]	1001
[5144]	1002	xenv1 = -(senv + deltasenv) / (sqrt(2.) * sigma1s)
	1003	xenv2 = -(senv - deltasenv) / (sqrt(2.) * sigma1s)
	1004	exp_xenv1 = exp(-1. * xenv1**2)
	1005	exp_xenv2 = exp(-1. * xenv2**2)
	1006	xth1 = -(sth + deltasth) / (sqrt(2.) * sigma2s)
	1007	xth2 = -(sth - deltasth) / (sqrt(2.) * sigma2s)
	1008	exp_xth1 = exp(-1. * xth1**2)
	1009	exp_xth2 = exp(-1. * xth2**2)
[2686]	1010
[5144]	1011	!CF_surfacique
	1012	cth(ind1, ind2) = 0.5 * (1. - 1. * erf(xth1))
	1013	cenv(ind1, ind2) = 0.5 * (1. - 1. * erf(xenv1))
	1014	ctot(ind1, ind2) = fraca(ind1, ind2) * cth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv(ind1, ind2)
[3493]	1015
	1016
[5144]	1017	!CF_volumique & eau condense
	1018	!environnement
	1019	IntJ = 0.5 * senv * (1 - erf(xenv2)) + (sigma1s / sqrt2pi) * exp_xenv2
	1020	IntJ_CF = 0.5 * (1. - 1. * erf(xenv2))
	1021	IF (deltasenv < 1.e-10) THEN
	1022	qlenv(ind1, ind2) = IntJ
	1023	cenv_vol(ind1, ind2) = IntJ_CF
	1024	else
	1025	IntI1 = (((senv + deltasenv)2 + (sigma1s)2) / (8 * deltasenv)) * (erf(xenv2) - erf(xenv1))
	1026	IntI2 = (sigma1s*2 / (4 deltasenv * sqrtpi)) * (xenv1 * exp_xenv1 - xenv2 * exp_xenv2)
	1027	IntI3 = ((sqrt2 * sigma1s * (senv + deltasenv)) / (4 * sqrtpi * deltasenv)) * (exp_xenv1 - exp_xenv2)
	1028	IntI1_CF = ((senv + deltasenv) * (erf(xenv2) - erf(xenv1))) / (4 * deltasenv)
	1029	IntI3_CF = (sqrt2 * sigma1s * (exp_xenv1 - exp_xenv2)) / (4 * sqrtpi * deltasenv)
	1030	qlenv(ind1, ind2) = IntJ + IntI1 + IntI2 + IntI3
	1031	cenv_vol(ind1, ind2) = IntJ_CF + IntI1_CF + IntI3_CF
	1032	endif
[3493]	1033
[5144]	1034	!thermique
	1035	IntJ = 0.5 * sth * (1 - erf(xth2)) + (sigma2s / sqrt2pi) * exp_xth2
	1036	IntJ_CF = 0.5 * (1. - 1. * erf(xth2))
	1037	IF (deltasth < 1.e-10) THEN
	1038	qlth(ind1, ind2) = IntJ
	1039	cth_vol(ind1, ind2) = IntJ_CF
	1040	else
	1041	IntI1 = (((sth + deltasth)2 + (sigma2s)2) / (8 * deltasth)) * (erf(xth2) - erf(xth1))
	1042	IntI2 = (sigma2s*2 / (4 deltasth * sqrtpi)) * (xth1 * exp_xth1 - xth2 * exp_xth2)
	1043	IntI3 = ((sqrt2 * sigma2s * (sth + deltasth)) / (4 * sqrtpi * deltasth)) * (exp_xth1 - exp_xth2)
	1044	IntI1_CF = ((sth + deltasth) * (erf(xth2) - erf(xth1))) / (4 * deltasth)
	1045	IntI3_CF = (sqrt2 * sigma2s * (exp_xth1 - exp_xth2)) / (4 * sqrtpi * deltasth)
	1046	qlth(ind1, ind2) = IntJ + IntI1 + IntI2 + IntI3
	1047	cth_vol(ind1, ind2) = IntJ_CF + IntI1_CF + IntI3_CF
	1048	endif
[3493]	1049
[5144]	1050	qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2)
	1051	ctot_vol(ind1, ind2) = fraca(ind1, ind2) * cth_vol(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv_vol(ind1, ind2)
[3493]	1052
[5144]	1053	ELSE IF (iflag_cloudth_vert == 5) THEN
	1054	sigma1s = (0.71794 + 0.000498239 * dz(ind1, ind2)) * (fraca(ind1, ind2)**0.5) &
	1055	/ (1 - fraca(ind1, ind2)) * (((sth - senv)2)0.5) &
	1056	+ ratqs(ind1, ind2) * po(ind1) !Environment
	1057	sigma2s = (0.03218 + 0.000092655 * dz(ind1, ind2)) / ((fraca(ind1, ind2) + 0.02)*0.5) (((sth - senv)2)0.5) + 0.002 * zqta(ind1, ind2) !Thermals
	1058	!sigma1s=(1.1*0.5)(fraca(ind1,ind2)*0.6)/(1-fraca(ind1,ind2))((sth-senv)2)0.5+0.002*po(ind1)
	1059	!sigma2s=0.11((sth-senv)2)0.5/(fraca(ind1,ind2)+0.01)0.4+0.002zqta(ind1,ind2)
	1060	xth = sth / (sqrt(2.) * sigma2s)
	1061	xenv = senv / (sqrt(2.) * sigma1s)
[3493]	1062
[5144]	1063	!Volumique
	1064	cth_vol(ind1, ind2) = 0.5 * (1. + 1. * erf(xth))
	1065	cenv_vol(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv))
	1066	ctot_vol(ind1, ind2) = fraca(ind1, ind2) * cth_vol(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv_vol(ind1, ind2)
[5160]	1067	!PRINT *,'jeanjean_CV=',ctot_vol(ind1,ind2)
[3493]	1068
[5144]	1069	qlth(ind1, ind2) = sigma2s * ((exp(-1. * xth*2) / sqrt2pi) + xth sqrt(2.) * cth_vol(ind1, ind2))
	1070	qlenv(ind1, ind2) = sigma1s * ((exp(-1. * xenv*2) / sqrt2pi) + xenv sqrt(2.) * cenv_vol(ind1, ind2))
	1071	qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2)
[3493]	1072
[5144]	1073	!Surfacique
	1074	!Neggers
	1075	!beta=0.0044
	1076	!inverse_rho=1.+beta*dz(ind1,ind2)
[5160]	1077	!PRINT *,'jeanjean : beta=',beta
[5144]	1078	!cth(ind1,ind2)=cth_vol(ind1,ind2)*inverse_rho
	1079	!cenv(ind1,ind2)=cenv_vol(ind1,ind2)*inverse_rho
	1080	!ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
[3493]	1081
[5144]	1082	!Brooks
	1083	a_Brooks = 0.6694
	1084	b_Brooks = 0.1882
	1085	A_Maj_Brooks = 0.1635 !-- sans shear
	1086	!A_Maj_Brooks=0.17 !-- ARM LES
	1087	!A_Maj_Brooks=0.18 !-- RICO LES
	1088	!A_Maj_Brooks=0.19 !-- BOMEX LES
	1089	Dx_Brooks = 200000.
	1090	f_Brooks = A_Maj_Brooks * (dz(ind1, ind2)*(a_Brooks)) (Dx_Brooks**(-b_Brooks))
[5160]	1091	!PRINT *,'jeanjean_f=',f_Brooks
[3493]	1092
[5144]	1093	cth(ind1, ind2) = 1. / (1. + exp(-1. * f_Brooks) * ((1. / max(1.e-15, min(cth_vol(ind1, ind2), 1.))) - 1.))
	1094	cenv(ind1, ind2) = 1. / (1. + exp(-1. * f_Brooks) * ((1. / max(1.e-15, min(cenv_vol(ind1, ind2), 1.))) - 1.))
	1095	ctot(ind1, ind2) = 1. / (1. + exp(-1. * f_Brooks) * ((1. / max(1.e-15, min(ctot_vol(ind1, ind2), 1.))) - 1.))
[5160]	1096	!PRINT *,'JJ_ctot_1',ctot(ind1,ind2)
[2686]	1097
[5144]	1098	ENDIF ! of if (iflag_cloudth_vert<5)
	1099	ENDIF ! of if (iflag_cloudth_vert==1 or 3 or 4)
[2686]	1100
[5144]	1101	! if (ctot(ind1,ind2).lt.1.e-10) THEN
	1102	IF (cenv(ind1, ind2)<1.e-10.OR.cth(ind1, ind2)<1.e-10) THEN
	1103	ctot(ind1, ind2) = 0.
	1104	ctot_vol(ind1, ind2) = 0.
	1105	qcloud(ind1) = zqsatenv(ind1, ind2)
[2686]	1106
[5144]	1107	else
[2686]	1108
[5144]	1109	qcloud(ind1) = qltot(ind1, ind2) / ctot(ind1, ind2) + zqs(ind1)
	1110	! qcloud(ind1)=fraca(ind1,ind2)*qlth(ind1,ind2)/cth(ind1,ind2) &
	1111	! & +(1.-1.fraca(ind1,ind2))qlenv(ind1,ind2)/cenv(ind1,ind2)+zqs(ind1)
	1112
	1113	endif
	1114
[3493]	1115	else ! gaussienne environnement seule
[3999]	1116
[5144]	1117	zqenv(ind1) = po(ind1)
	1118	Tbef = t(ind1, ind2)
	1119	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	1120	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	1121	qsatbef = MIN(0.5, qsatbef)
	1122	zcor = 1. / (1. - retv * qsatbef)
	1123	qsatbef = qsatbef * zcor
	1124	zqsatenv(ind1, ind2) = qsatbef
[2686]	1125
	1126
[5144]	1127	! qlbef=Max(po(ind1)-zqsatenv(ind1,ind2),0.)
	1128	zthl(ind1, ind2) = t(ind1, ind2) * (101325 / paprs(ind1, ind2))**(rdd / cppd)
	1129	alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2)
	1130	aenv = 1. / (1. + (alenv * Lv / cppd))
	1131	senv = aenv * (po(ind1) - zqsatenv(ind1, ind2))
	1132	sth = 0.
[2686]	1133
[5144]	1134	sigma1s = ratqs(ind1, ind2) * zqenv(ind1)
	1135	sigma2s = 0.
[2686]	1136
[5144]	1137	sqrt2pi = sqrt(2. * pi)
	1138	xenv = senv / (sqrt(2.) * sigma1s)
	1139	ctot(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv))
	1140	ctot_vol(ind1, ind2) = ctot(ind1, ind2)
	1141	qltot(ind1, ind2) = sigma1s * ((exp(-1. * xenv*2) / sqrt2pi) + xenv sqrt(2.) * cenv(ind1, ind2))
[2686]	1142
[5144]	1143	IF (ctot(ind1, ind2)<1.e-3) THEN
	1144	ctot(ind1, ind2) = 0.
	1145	qcloud(ind1) = zqsatenv(ind1, ind2)
[3493]	1146
[5144]	1147	else
[3493]	1148
[5144]	1149	! ctot(ind1,ind2)=ctot(ind1,ind2)
	1150	qcloud(ind1) = qltot(ind1, ind2) / ctot(ind1, ind2) + zqsatenv(ind1, ind2)
[3493]	1151
[5144]	1152	endif
	1153
[2686]	1154	endif ! From the separation (thermal/envrionnement) et (environnement) only, l.335 et l.492
[2958]	1155	! Outputs used to check the PDFs
[5144]	1156	cloudth_senv(ind1, ind2) = senv
	1157	cloudth_sth(ind1, ind2) = sth
	1158	cloudth_sigmaenv(ind1, ind2) = sigma1s
	1159	cloudth_sigmath(ind1, ind2) = sigma2s
[2958]	1160
[5144]	1161	enddo ! from the loop on ngrid l.333
	1162	RETURN
	1163	! end
	1164	END SUBROUTINE cloudth_vert_v3
[3493]	1165
[5144]	1166	SUBROUTINE cloudth_v6(ngrid, klev, ind2, &
	1167	ztv, po, zqta, fraca, &
	1168	qcloud, ctot_surf, ctot_vol, zpspsk, paprs, pplay, ztla, zthl, &
	1169	ratqs, zqs, T, &
	1170	cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv)
[3493]	1171
[5144]	1172	USE lmdz_cloudth_ini, ONLY: iflag_cloudth_vert
	1173	USE lmdz_yoethf
[5153]	1174
[5144]	1175	USE lmdz_yomcst
[3493]	1176
[5144]	1177	IMPLICIT NONE
[5153]	1178	INCLUDE "FCTTRE.h"
[3493]	1179
[5144]	1180	!Domain variables
	1181	INTEGER ngrid !indice Max lat-lon
	1182	INTEGER klev !indice Max alt
	1183	REAL, DIMENSION(ngrid, klev), INTENT(OUT) :: cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv
	1184	INTEGER ind1 !indice in [1:ngrid]
	1185	INTEGER ind2 !indice in [1:klev]
	1186	!thermal plume fraction
	1187	REAL fraca(ngrid, klev + 1) !thermal plumes fraction in the gridbox
	1188	!temperatures
	1189	REAL T(ngrid, klev) !temperature
	1190	REAL zpspsk(ngrid, klev) !factor (p/p0)**kappa (used for potential variables)
	1191	REAL ztv(ngrid, klev) !potential temperature (voir thermcell_env.F90)
	1192	REAL ztla(ngrid, klev) !liquid temperature in the thermals (Tl_th)
	1193	REAL zthl(ngrid, klev) !liquid temperature in the environment (Tl_env)
	1194	!pressure
	1195	REAL paprs(ngrid, klev + 1) !pressure at the interface of levels
	1196	REAL pplay(ngrid, klev) !pressure at the middle of the level
	1197	!humidity
	1198	REAL ratqs(ngrid, klev) !width of the total water subgrid-scale distribution
	1199	REAL po(ngrid) !total water (qt)
	1200	REAL zqenv(ngrid) !total water in the environment (qt_env)
	1201	REAL zqta(ngrid, klev) !total water in the thermals (qt_th)
	1202	REAL zqsatth(ngrid, klev) !water saturation level in the thermals (q_sat_th)
	1203	REAL zqsatenv(ngrid, klev) !water saturation level in the environment (q_sat_env)
	1204	REAL qlth(ngrid, klev) !condensed water in the thermals
	1205	REAL qlenv(ngrid, klev) !condensed water in the environment
	1206	REAL qltot(ngrid, klev) !condensed water in the gridbox
	1207	!cloud fractions
	1208	REAL cth_vol(ngrid, klev) !cloud fraction by volume in the thermals
	1209	REAL cenv_vol(ngrid, klev) !cloud fraction by volume in the environment
	1210	REAL ctot_vol(ngrid, klev) !cloud fraction by volume in the gridbox
	1211	REAL cth_surf(ngrid, klev) !cloud fraction by surface in the thermals
	1212	REAL cenv_surf(ngrid, klev) !cloud fraction by surface in the environment
	1213	REAL ctot_surf(ngrid, klev) !cloud fraction by surface in the gridbox
	1214	!PDF of saturation deficit variables
	1215	REAL rdd, cppd, Lv
	1216	REAL Tbef, zdelta, qsatbef, zcor
	1217	REAL alth, alenv, ath, aenv
	1218	REAL sth, senv !saturation deficits in the thermals and environment
	1219	REAL sigma_env, sigma_th !standard deviations of the biGaussian PDF
	1220	!cloud fraction variables
	1221	REAL xth, xenv
	1222	REAL inverse_rho, beta !Neggers et al. (2011) method
	1223	REAL a_Brooks, b_Brooks, A_Maj_Brooks, Dx_Brooks, f_Brooks !Brooks et al. (2005) method
	1224	!Incloud total water variables
	1225	REAL zqs(ngrid) !q_sat
	1226	REAL qcloud(ngrid) !eau totale dans le nuage
	1227	!Some arithmetic variables
	1228	REAL erf, pi, sqrt2, sqrt2pi
	1229	!Depth of the layer
	1230	REAL dz(ngrid, klev) !epaisseur de la couche en metre
	1231	REAL rhodz(ngrid, klev)
	1232	REAL zrho(ngrid, klev)
	1233	DO ind1 = 1, ngrid
	1234	rhodz(ind1, ind2) = (paprs(ind1, ind2) - paprs(ind1, ind2 + 1)) / rg ![kg/m2]
	1235	zrho(ind1, ind2) = pplay(ind1, ind2) / T(ind1, ind2) / rd ![kg/m3]
	1236	dz(ind1, ind2) = rhodz(ind1, ind2) / zrho(ind1, ind2) ![m]
	1237	END DO
[3493]	1238
[5144]	1239	!------------------------------------------------------------------------------
	1240	! Initialization
	1241	!------------------------------------------------------------------------------
	1242	qlth(:, :) = 0.
	1243	qlenv(:, :) = 0.
	1244	qltot(:, :) = 0.
	1245	cth_vol(:, :) = 0.
	1246	cenv_vol(:, :) = 0.
	1247	ctot_vol(:, :) = 0.
	1248	cth_surf(:, :) = 0.
	1249	cenv_surf(:, :) = 0.
	1250	ctot_surf(:, :) = 0.
	1251	qcloud(:) = 0.
	1252	rdd = 287.04
	1253	cppd = 1005.7
	1254	pi = 3.14159
	1255	Lv = 2.5e6
	1256	sqrt2 = sqrt(2.)
	1257	sqrt2pi = sqrt(2. * pi)
[3493]	1258
[5144]	1259	DO ind1 = 1, ngrid
	1260	!-------------------------------------------------------------------------------
	1261	!Both thermal and environment in the gridbox
	1262	!-------------------------------------------------------------------------------
	1263	IF ((ztv(ind1, 1)>ztv(ind1, 2)).AND.(fraca(ind1, ind2)>1.e-10)) THEN
[3493]	1264	!--------------------------------------------
	1265	!calcul de qsat_env
	1266	!--------------------------------------------
[5144]	1267	Tbef = zthl(ind1, ind2) * zpspsk(ind1, ind2)
	1268	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	1269	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	1270	qsatbef = MIN(0.5, qsatbef)
	1271	zcor = 1. / (1. - retv * qsatbef)
	1272	qsatbef = qsatbef * zcor
	1273	zqsatenv(ind1, ind2) = qsatbef
[3493]	1274	!--------------------------------------------
	1275	!calcul de s_env
	1276	!--------------------------------------------
[5144]	1277	alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2) !qsl, p84 these Arnaud Jam
	1278	aenv = 1. / (1. + (alenv * Lv / cppd)) !al, p84 these Arnaud Jam
	1279	senv = aenv * (po(ind1) - zqsatenv(ind1, ind2)) !s, p84 these Arnaud Jam
[3493]	1280	!--------------------------------------------
	1281	!calcul de qsat_th
	1282	!--------------------------------------------
[5144]	1283	Tbef = ztla(ind1, ind2) * zpspsk(ind1, ind2)
	1284	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	1285	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	1286	qsatbef = MIN(0.5, qsatbef)
	1287	zcor = 1. / (1. - retv * qsatbef)
	1288	qsatbef = qsatbef * zcor
	1289	zqsatth(ind1, ind2) = qsatbef
[3493]	1290	!--------------------------------------------
[5144]	1291	!calcul de s_th
[3493]	1292	!--------------------------------------------
[5144]	1293	alth = (0.622 * Lv * zqsatth(ind1, ind2)) / (rdd * ztla(ind1, ind2)**2) !qsl, p84 these Arnaud Jam
	1294	ath = 1. / (1. + (alth * Lv / cppd)) !al, p84 these Arnaud Jam
	1295	sth = ath * (zqta(ind1, ind2) - zqsatth(ind1, ind2)) !s, p84 these Arnaud Jam
[3493]	1296	!--------------------------------------------
	1297	!calcul standard deviations bi-Gaussian PDF
	1298	!--------------------------------------------
[5144]	1299	sigma_th = (0.03218 + 0.000092655 * dz(ind1, ind2)) / ((fraca(ind1, ind2) + 0.01)*0.5) (((sth - senv)2)0.5) + 0.002 * zqta(ind1, ind2)
	1300	sigma_env = (0.71794 + 0.000498239 * dz(ind1, ind2)) * (fraca(ind1, ind2)**0.5) &
	1301	/ (1 - fraca(ind1, ind2)) * (((sth - senv)2)0.5) &
	1302	+ ratqs(ind1, ind2) * po(ind1)
	1303	xth = sth / (sqrt2 * sigma_th)
	1304	xenv = senv / (sqrt2 * sigma_env)
[3493]	1305	!--------------------------------------------
	1306	!Cloud fraction by volume CF_vol
	1307	!--------------------------------------------
[5144]	1308	cth_vol(ind1, ind2) = 0.5 * (1. + 1. * erf(xth))
	1309	cenv_vol(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv))
	1310	ctot_vol(ind1, ind2) = fraca(ind1, ind2) * cth_vol(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * cenv_vol(ind1, ind2)
	1311	!--------------------------------------------
[3493]	1312	!Condensed water qc
	1313	!--------------------------------------------
[5144]	1314	qlth(ind1, ind2) = sigma_th * ((exp(-1. * xth*2) / sqrt2pi) + xth sqrt2 * cth_vol(ind1, ind2))
	1315	qlenv(ind1, ind2) = sigma_env * ((exp(-1. * xenv*2) / sqrt2pi) + xenv sqrt2 * cenv_vol(ind1, ind2))
	1316	qltot(ind1, ind2) = fraca(ind1, ind2) * qlth(ind1, ind2) + (1. - 1. * fraca(ind1, ind2)) * qlenv(ind1, ind2)
[3493]	1317	!--------------------------------------------
	1318	!Cloud fraction by surface CF_surf
	1319	!--------------------------------------------
	1320	!Method Neggers et al. (2011) : ok for cumulus clouds only
[5144]	1321	!beta=0.0044 (Jouhaud et al.2018)
	1322	!inverse_rho=1.+beta*dz(ind1,ind2)
	1323	!ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
[3493]	1324	!Method Brooks et al. (2005) : ok for all types of clouds
[5144]	1325	a_Brooks = 0.6694
	1326	b_Brooks = 0.1882
	1327	A_Maj_Brooks = 0.1635 !-- sans dependence au cisaillement de vent
	1328	Dx_Brooks = 200000. !-- si l'on considere des mailles de 200km de cote
	1329	f_Brooks = A_Maj_Brooks * (dz(ind1, ind2)*(a_Brooks)) (Dx_Brooks**(-b_Brooks))
	1330	ctot_surf(ind1, ind2) = 1. / (1. + exp(-1. * f_Brooks) * ((1. / max(1.e-15, min(ctot_vol(ind1, ind2), 1.))) - 1.))
[3493]	1331	!--------------------------------------------
	1332	!Incloud Condensed water qcloud
	1333	!--------------------------------------------
[5144]	1334	IF (ctot_surf(ind1, ind2) < 1.e-10) THEN
	1335	ctot_vol(ind1, ind2) = 0.
	1336	ctot_surf(ind1, ind2) = 0.
	1337	qcloud(ind1) = zqsatenv(ind1, ind2)
	1338	else
	1339	qcloud(ind1) = qltot(ind1, ind2) / ctot_vol(ind1, ind2) + zqs(ind1)
	1340	endif
[3493]	1341
	1342
	1343
[5144]	1344	!-------------------------------------------------------------------------------
	1345	!Environment only in the gridbox
	1346	!-------------------------------------------------------------------------------
[3493]	1347	ELSE
	1348	!--------------------------------------------
	1349	!calcul de qsat_env
	1350	!--------------------------------------------
[5144]	1351	Tbef = zthl(ind1, ind2) * zpspsk(ind1, ind2)
	1352	zdelta = MAX(0., SIGN(1., RTT - Tbef))
	1353	qsatbef = R2ES * FOEEW(Tbef, zdelta) / paprs(ind1, ind2)
	1354	qsatbef = MIN(0.5, qsatbef)
	1355	zcor = 1. / (1. - retv * qsatbef)
	1356	qsatbef = qsatbef * zcor
	1357	zqsatenv(ind1, ind2) = qsatbef
[3493]	1358	!--------------------------------------------
	1359	!calcul de s_env
	1360	!--------------------------------------------
[5144]	1361	alenv = (0.622 * Lv * zqsatenv(ind1, ind2)) / (rdd * zthl(ind1, ind2)**2) !qsl, p84 these Arnaud Jam
	1362	aenv = 1. / (1. + (alenv * Lv / cppd)) !al, p84 these Arnaud Jam
	1363	senv = aenv * (po(ind1) - zqsatenv(ind1, ind2)) !s, p84 these Arnaud Jam
[3493]	1364	!--------------------------------------------
	1365	!calcul standard deviations Gaussian PDF
	1366	!--------------------------------------------
[5144]	1367	zqenv(ind1) = po(ind1)
	1368	sigma_env = ratqs(ind1, ind2) * zqenv(ind1)
	1369	xenv = senv / (sqrt2 * sigma_env)
[3493]	1370	!--------------------------------------------
	1371	!Cloud fraction by volume CF_vol
	1372	!--------------------------------------------
[5144]	1373	ctot_vol(ind1, ind2) = 0.5 * (1. + 1. * erf(xenv))
[3493]	1374	!--------------------------------------------
	1375	!Condensed water qc
	1376	!--------------------------------------------
[5144]	1377	qltot(ind1, ind2) = sigma_env * ((exp(-1. * xenv*2) / sqrt2pi) + xenv sqrt2 * ctot_vol(ind1, ind2))
[3493]	1378	!--------------------------------------------
	1379	!Cloud fraction by surface CF_surf
	1380	!--------------------------------------------
	1381	!Method Neggers et al. (2011) : ok for cumulus clouds only
[5144]	1382	!beta=0.0044 (Jouhaud et al.2018)
	1383	!inverse_rho=1.+beta*dz(ind1,ind2)
	1384	!ctot_surf(ind1,ind2)=ctot_vol(ind1,ind2)*inverse_rho
[3493]	1385	!Method Brooks et al. (2005) : ok for all types of clouds
[5144]	1386	a_Brooks = 0.6694
	1387	b_Brooks = 0.1882
	1388	A_Maj_Brooks = 0.1635 !-- sans dependence au shear
	1389	Dx_Brooks = 200000.
	1390	f_Brooks = A_Maj_Brooks * (dz(ind1, ind2)*(a_Brooks)) (Dx_Brooks**(-b_Brooks))
	1391	ctot_surf(ind1, ind2) = 1. / (1. + exp(-1. * f_Brooks) * ((1. / max(1.e-15, min(ctot_vol(ind1, ind2), 1.))) - 1.))
[3493]	1392	!--------------------------------------------
	1393	!Incloud Condensed water qcloud
	1394	!--------------------------------------------
[5144]	1395	IF (ctot_surf(ind1, ind2) < 1.e-8) THEN
	1396	ctot_vol(ind1, ind2) = 0.
	1397	ctot_surf(ind1, ind2) = 0.
	1398	qcloud(ind1) = zqsatenv(ind1, ind2)
	1399	else
	1400	qcloud(ind1) = qltot(ind1, ind2) / ctot_vol(ind1, ind2) + zqsatenv(ind1, ind2)
	1401	endif
[3493]	1402
	1403	END IF ! From the separation (thermal/envrionnement) et (environnement only)
	1404
	1405	! Outputs used to check the PDFs
[5144]	1406	cloudth_senv(ind1, ind2) = senv
	1407	cloudth_sth(ind1, ind2) = sth
	1408	cloudth_sigmaenv(ind1, ind2) = sigma_env
	1409	cloudth_sigmath(ind1, ind2) = sigma_th
[3493]	1410
[5144]	1411	END DO ! From the loop on ngrid
[3493]	1412
[5144]	1413	END SUBROUTINE cloudth_v6
[5105]	1414
[3999]	1415
[5144]	1416	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	1417	SUBROUTINE cloudth_mpc(klon, klev, ind2, mpc_bl_points, &
	1418	& temp, qt, qt_th, frac_th, zpspsk, paprsup, paprsdn, play, thetal_th, &
	1419	& ratqs, qcloud, qincloud, icefrac, ctot, ctot_vol, &
	1420	& cloudth_sth, cloudth_senv, cloudth_sigmath, cloudth_sigmaenv)
	1421	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	1422	! Author : Etienne Vignon (LMDZ/CNRS)
	1423	! Date: April 2024
	1424	! Date: Adapted from cloudth_vert_v3 in 2023 by Arnaud Otavio Jam
	1425	! IMPORTANT NOTE: we assume iflag_cloudth_vert=7
	1426	!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
[3999]	1427
[5144]	1428	USE lmdz_cloudth_ini, ONLY: iflag_cloudth_vert, iflag_ratqs
	1429	USE lmdz_cloudth_ini, ONLY: C_mpc, Ni, C_cap, Ei, d_top, vert_alpha, vert_alpha_th, sigma1s_factor, sigma1s_power, sigma2s_factor, sigma2s_power, cloudth_ratqsmin, iflag_cloudth_vert_noratqs
	1430	USE lmdz_lscp_tools, only: CALC_QSAT_ECMWF
	1431	USE lmdz_lscp_ini, only: RTT, RG, RPI, RD, RCPD, RLVTT, RLSTT, temp_nowater, min_frac_th_cld, min_neb_th
[3999]	1432
[5144]	1433	IMPLICIT NONE
[3999]	1434
	1435
[5144]	1436	!------------------------------------------------------------------------------
	1437	! Declaration
	1438	!------------------------------------------------------------------------------
[3999]	1439
[5144]	1440	! INPUT/OUTPUT
[3999]	1441
[5144]	1442	INTEGER, INTENT(IN) :: klon, klev, ind2
[3999]	1443
[5144]	1444	REAL, DIMENSION(klon), INTENT(IN) :: temp ! Temperature (liquid temperature) in the mesh [K] : has seen evap of precip
	1445	REAL, DIMENSION(klon), INTENT(IN) :: qt ! total water specific humidity in the mesh [kg/kg]: has seen evap of precip
	1446	REAL, DIMENSION(klon), INTENT(IN) :: qt_th ! total water specific humidity in thermals [kg/kg]: has not seen evap of precip
	1447	REAL, DIMENSION(klon), INTENT(IN) :: thetal_th ! Liquid potential temperature in thermals [K]: has not seen the evap of precip
	1448	REAL, DIMENSION(klon), INTENT(IN) :: frac_th ! Fraction of the mesh covered by thermals [0-1]
	1449	REAL, DIMENSION(klon), INTENT(IN) :: zpspsk ! Exner potential
	1450	REAL, DIMENSION(klon), INTENT(IN) :: paprsup ! Pressure at top layer interface [Pa]
	1451	REAL, DIMENSION(klon), INTENT(IN) :: paprsdn ! Pressure at bottom layer interface [Pa]
	1452	REAL, DIMENSION(klon), INTENT(IN) :: play ! Pressure of layers [Pa]
	1453	REAL, DIMENSION(klon), INTENT(IN) :: ratqs ! Parameter that determines the width of the total water distrib.
[3999]	1454
[5144]	1455	INTEGER, DIMENSION(klon, klev), INTENT(INOUT) :: mpc_bl_points ! grid points with convective (thermals) mixed phase clouds
[3999]	1456
[5144]	1457	REAL, DIMENSION(klon), INTENT(OUT) :: ctot ! Cloud fraction [0-1]
	1458	REAL, DIMENSION(klon), INTENT(OUT) :: ctot_vol ! Volume cloud fraction [0-1]
	1459	REAL, DIMENSION(klon), INTENT(OUT) :: qcloud ! In cloud total water content [kg/kg]
	1460	REAL, DIMENSION(klon), INTENT(OUT) :: qincloud ! In cloud condensed water content [kg/kg]
	1461	REAL, DIMENSION(klon), INTENT(OUT) :: icefrac ! Fraction of ice in clouds [0-1]
	1462	REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_sth ! mean saturation deficit in thermals
	1463	REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_senv ! mean saturation deficit in environment
	1464	REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_sigmath ! std of saturation deficit in thermals
	1465	REAL, DIMENSION(klon), INTENT(OUT) :: cloudth_sigmaenv ! std of saturation deficit in environment
[3999]	1466
	1467
[5144]	1468	! LOCAL VARIABLES
[3999]	1469
[5144]	1470	INTEGER itap, ind1, l, ig, iter, k
	1471	INTEGER iflag_topthermals, niter
[4380]	1472
[5144]	1473	REAL qcth(klon)
	1474	REAL qcenv(klon)
	1475	REAL qctot(klon)
	1476	REAL cth(klon)
	1477	REAL cenv(klon)
	1478	REAL cth_vol(klon)
	1479	REAL cenv_vol(klon)
	1480	REAL qt_env(klon), thetal_env(klon)
	1481	REAL icefracenv(klon), icefracth(klon)
	1482	REAL sqrtpi, sqrt2, sqrt2pi
	1483	REAL alth, alenv, ath, aenv
	1484	REAL sth, senv, sigma1s, sigma2s, sigma1s_fraca, sigma1s_ratqs
	1485	REAL inverse_rho, beta, a_Brooks, b_Brooks, A_Maj_Brooks, Dx_Brooks, f_Brooks
	1486	REAL xth, xenv, exp_xenv1, exp_xenv2, exp_xth1, exp_xth2
	1487	REAL xth1, xth2, xenv1, xenv2, deltasth, deltasenv
	1488	REAL IntJ, IntI1, IntI2, IntI3, IntJ_CF, IntI1_CF, IntI3_CF, coeffqlenv, coeffqlth
	1489	REAL zdelta, qsatbef, zcor
	1490	REAL Tbefth(klon), Tbefenv(klon), zrho(klon), rhoth(klon)
	1491	REAL qlbef
	1492	REAL dqsatenv(klon), dqsatth(klon)
	1493	REAL erf
	1494	REAL zpdf_sig(klon), zpdf_k(klon), zpdf_delta(klon)
	1495	REAL zpdf_a(klon), zpdf_b(klon), zpdf_e1(klon), zpdf_e2(klon)
	1496	REAL rhodz(klon)
	1497	REAL rho(klon)
	1498	REAL dz(klon)
	1499	REAL qslenv(klon), qslth(klon), qsienv(klon), qsith(klon)
	1500	REAL alenvl, aenvl
	1501	REAL sthi, sthl, sthil, althl, athl, althi, athi, sthlc, deltasthc, sigma2sc
	1502	REAL senvi, senvl, qbase, sbase, qliqth, qiceth
	1503	REAL qmax, ttarget, stmp, cout, coutref, fraci
	1504	REAL maxi, mini, pas
[3999]	1505
[5144]	1506	! Modifty the saturation deficit PDF in thermals
	1507	! in the presence of ice crystals
	1508	CHARACTER (len = 80) :: abort_message
	1509	CHARACTER (len = 20) :: routname = 'cloudth_mpc'
[3999]	1510
	1511
[5144]	1512	!------------------------------------------------------------------------------
	1513	! Initialisation
	1514	!------------------------------------------------------------------------------
[3999]	1515
	1516
[5144]	1517	! Few initial checksS
[3999]	1518
[5144]	1519	DO k = 1, klev
[3999]	1520	DO ind1 = 1, klon
[5144]	1521	rhodz(ind1) = (paprsdn(ind1) - paprsup(ind1)) / rg !kg/m2
	1522	zrho(ind1) = play(ind1) / temp(ind1) / rd !kg/m3
	1523	dz(ind1) = rhodz(ind1) / zrho(ind1) !m : epaisseur de la couche en metre
[3999]	1524	END DO
[5144]	1525	END DO
[3999]	1526
[5144]	1527	icefracth(:) = 0.
	1528	icefracenv(:) = 0.
	1529	sqrt2pi = sqrt(2. * rpi)
	1530	sqrt2 = sqrt(2.)
	1531	sqrtpi = sqrt(rpi)
	1532	icefrac(:) = 0.
[3999]	1533
[5144]	1534	!-------------------------------------------------------------------------------
	1535	! Identify grid points with potential mixed-phase conditions
	1536	!-------------------------------------------------------------------------------
[3999]	1537
[5144]	1538	DO ind1 = 1, klon
	1539	IF ((temp(ind1) < RTT) .AND. (temp(ind1) > temp_nowater) &
	1540	.AND. (ind2<=klev - 2) &
	1541	.AND. (frac_th(ind1)>min_frac_th_cld)) THEN
	1542	mpc_bl_points(ind1, ind2) = 1
	1543	ELSE
	1544	mpc_bl_points(ind1, ind2) = 0
	1545	ENDIF
	1546	ENDDO
[3999]	1547
	1548
[5144]	1549	!-------------------------------------------------------------------------------
	1550	! Thermal fraction calculation and standard deviation of the distribution
	1551	!-------------------------------------------------------------------------------
[3999]	1552
[5144]	1553	! initialisations and calculation of temperature, humidity and saturation specific humidity
[3999]	1554
[5144]	1555	DO ind1 = 1, klon
[3999]	1556
[5144]	1557	rhodz(ind1) = (paprsdn(ind1) - paprsup(ind1)) / rg ! kg/m2
	1558	rho(ind1) = play(ind1) / temp(ind1) / rd ! kg/m3
	1559	dz(ind1) = rhodz(ind1) / rho(ind1) ! m : epaisseur de la couche en metre
	1560	Tbefenv(ind1) = temp(ind1)
	1561	thetal_env(ind1) = Tbefenv(ind1) / zpspsk(ind1)
	1562	Tbefth(ind1) = thetal_th(ind1) * zpspsk(ind1)
	1563	rhoth(ind1) = play(ind1) / Tbefth(ind1) / RD
	1564	qt_env(ind1) = (qt(ind1) - frac_th(ind1) * qt_th(ind1)) / (1. - frac_th(ind1)) !qt = aqtth + (1-a)qtenv
[3999]	1565
[5144]	1566	ENDDO
[3999]	1567
[5144]	1568	! Calculation of saturation specific humidity
	1569	CALL CALC_QSAT_ECMWF(klon, Tbefenv, qt_env, play, RTT, 1, .FALSE., qslenv, dqsatenv)
	1570	CALL CALC_QSAT_ECMWF(klon, Tbefenv, qt_env, play, RTT, 2, .FALSE., qsienv, dqsatenv)
	1571	CALL CALC_QSAT_ECMWF(klon, Tbefth, qt_th, play, RTT, 1, .FALSE., qslth, dqsatth)
[4072]	1572
[5144]	1573	DO ind1 = 1, klon
[4072]	1574
[5144]	1575	IF (frac_th(ind1)>min_frac_th_cld) THEN !Thermal and environnement
[4072]	1576
[5144]	1577	! unlike in the other cloudth routine,
	1578	! We consider distributions of the saturation deficit WRT liquid
	1579	! at positive AND negative celcius temperature
	1580	! subsequently, cloud fraction corresponds to the part of the pdf corresponding
	1581	! to superstauration with respect to liquid WHATEVER the temperature
[4072]	1582
[5144]	1583	! Environment:
[4072]	1584
[5144]	1585	alenv = (0.622 * RLVTT * qslenv(ind1)) / (rd * thetal_env(ind1)**2)
	1586	aenv = 1. / (1. + (alenv * RLVTT / rcpd))
	1587	senv = aenv * (qt_env(ind1) - qslenv(ind1))
[4072]	1588
	1589
[5144]	1590	! Thermals:
[3999]	1591
[5144]	1592	alth = (0.622 * RLVTT * qslth(ind1)) / (rd * thetal_th(ind1)**2)
	1593	ath = 1. / (1. + (alth * RLVTT / rcpd))
	1594	sth = ath * (qt_th(ind1) - qslth(ind1))
[3999]	1595
	1596
[5144]	1597	! IF (mpc_bl_points(ind1,ind2) .EQ. 0) THEN ! No BL MPC
[3999]	1598
	1599
[5144]	1600	! Standard deviation of the distributions
[3999]	1601
[5144]	1602	sigma1s_fraca = (sigma1s_factor*0.5) (frac_th(ind1)**sigma1s_power) / &
	1603	(1 - frac_th(ind1)) * ((sth - senv)2)0.5
[3999]	1604
[5144]	1605	IF (cloudth_ratqsmin>0.) THEN
	1606	sigma1s_ratqs = cloudth_ratqsmin * qt(ind1)
	1607	ELSE
	1608	sigma1s_ratqs = ratqs(ind1) * qt(ind1)
	1609	ENDIF
[3999]	1610
[5144]	1611	sigma1s = sigma1s_fraca + sigma1s_ratqs
	1612	IF (iflag_ratqs==11) THEN
	1613	sigma1s = ratqs(ind1) * qt(ind1) * aenv
	1614	ENDIF
	1615	sigma2s = (sigma2s_factor * (((sth - senv)2)0.5) / ((frac_th(ind1) + 0.02)*sigma2s_power)) + 0.002 qt_th(ind1)
[3999]	1616
[5144]	1617	deltasenv = aenv * vert_alpha * sigma1s
	1618	deltasth = ath * vert_alpha_th * sigma2s
[3999]	1619
[5144]	1620	xenv1 = -(senv + deltasenv) / (sqrt(2.) * sigma1s)
	1621	xenv2 = -(senv - deltasenv) / (sqrt(2.) * sigma1s)
	1622	exp_xenv1 = exp(-1. * xenv1**2)
	1623	exp_xenv2 = exp(-1. * xenv2**2)
	1624	xth1 = -(sth + deltasth) / (sqrt(2.) * sigma2s)
	1625	xth2 = -(sth - deltasth) / (sqrt(2.) * sigma2s)
	1626	exp_xth1 = exp(-1. * xth1**2)
	1627	exp_xth2 = exp(-1. * xth2**2)
[3999]	1628
[5144]	1629	!surface CF
[3999]	1630
[5144]	1631	cth(ind1) = 0.5 * (1. - 1. * erf(xth1))
	1632	cenv(ind1) = 0.5 * (1. - 1. * erf(xenv1))
	1633	ctot(ind1) = frac_th(ind1) * cth(ind1) + (1. - 1. * frac_th(ind1)) * cenv(ind1)
[3999]	1634
	1635
[5144]	1636	!volume CF, condensed water and ice fraction
[3999]	1637
[5144]	1638	!environnement
[3999]	1639
[5144]	1640	IntJ = 0.5 * senv * (1 - erf(xenv2)) + (sigma1s / sqrt2pi) * exp_xenv2
	1641	IntJ_CF = 0.5 * (1. - 1. * erf(xenv2))
[3999]	1642
[5144]	1643	IF (deltasenv < 1.e-10) THEN
	1644	qcenv(ind1) = IntJ
	1645	cenv_vol(ind1) = IntJ_CF
	1646	ELSE
	1647	IntI1 = (((senv + deltasenv)2 + (sigma1s)2) / (8 * deltasenv)) * (erf(xenv2) - erf(xenv1))
	1648	IntI2 = (sigma1s*2 / (4 deltasenv * sqrtpi)) * (xenv1 * exp_xenv1 - xenv2 * exp_xenv2)
	1649	IntI3 = ((sqrt2 * sigma1s * (senv + deltasenv)) / (4 * sqrtpi * deltasenv)) * (exp_xenv1 - exp_xenv2)
	1650	IntI1_CF = ((senv + deltasenv) * (erf(xenv2) - erf(xenv1))) / (4 * deltasenv)
	1651	IntI3_CF = (sqrt2 * sigma1s * (exp_xenv1 - exp_xenv2)) / (4 * sqrtpi * deltasenv)
	1652	qcenv(ind1) = IntJ + IntI1 + IntI2 + IntI3
	1653	cenv_vol(ind1) = IntJ_CF + IntI1_CF + IntI3_CF
	1654	IF (Tbefenv(ind1) < temp_nowater) THEN
	1655	! freeze all droplets in cirrus temperature regime
	1656	icefracenv(ind1) = 1.
	1657	ENDIF
	1658	ENDIF
[3999]	1659
	1660
	1661
[5144]	1662	!thermals
[3999]	1663
[5144]	1664	IntJ = 0.5 * sth * (1 - erf(xth2)) + (sigma2s / sqrt2pi) * exp_xth2
	1665	IntJ_CF = 0.5 * (1. - 1. * erf(xth2))
[3999]	1666
[5144]	1667	IF (deltasth < 1.e-10) THEN
	1668	qcth(ind1) = IntJ
	1669	cth_vol(ind1) = IntJ_CF
	1670	ELSE
	1671	IntI1 = (((sth + deltasth)2 + (sigma2s)2) / (8 * deltasth)) * (erf(xth2) - erf(xth1))
	1672	IntI2 = (sigma2s*2 / (4 deltasth * sqrtpi)) * (xth1 * exp_xth1 - xth2 * exp_xth2)
	1673	IntI3 = ((sqrt2 * sigma2s * (sth + deltasth)) / (4 * sqrtpi * deltasth)) * (exp_xth1 - exp_xth2)
	1674	IntI1_CF = ((sth + deltasth) * (erf(xth2) - erf(xth1))) / (4 * deltasth)
	1675	IntI3_CF = (sqrt2 * sigma2s * (exp_xth1 - exp_xth2)) / (4 * sqrtpi * deltasth)
	1676	qcth(ind1) = IntJ + IntI1 + IntI2 + IntI3
	1677	cth_vol(ind1) = IntJ_CF + IntI1_CF + IntI3_CF
	1678	IF (Tbefth(ind1) < temp_nowater) THEN
	1679	! freeze all droplets in cirrus temperature regime
	1680	icefracth(ind1) = 1.
	1681	ENDIF
[3999]	1682
[5144]	1683	ENDIF
[4910]	1684
[5144]	1685	qctot(ind1) = frac_th(ind1) * qcth(ind1) + (1. - 1. * frac_th(ind1)) * qcenv(ind1)
	1686	ctot_vol(ind1) = frac_th(ind1) * cth_vol(ind1) + (1. - 1. * frac_th(ind1)) * cenv_vol(ind1)
[3999]	1687
[5144]	1688	IF (cenv(ind1)<min_neb_th.AND.cth(ind1)<min_neb_th) THEN
	1689	ctot(ind1) = 0.
	1690	ctot_vol(ind1) = 0.
	1691	qcloud(ind1) = qslenv(ind1)
	1692	qincloud(ind1) = 0.
	1693	icefrac(ind1) = 0.
	1694	ELSE
	1695	qcloud(ind1) = qctot(ind1) / ctot(ind1) + qslenv(ind1)
	1696	qincloud(ind1) = qctot(ind1) / ctot(ind1)
	1697	IF (qctot(ind1) > 0) THEN
	1698	icefrac(ind1) = (frac_th(ind1) * qcth(ind1) * icefracth(ind1) + (1. - 1. * frac_th(ind1)) * qcenv(ind1) * icefracenv(ind1)) / qctot(ind1)
	1699	icefrac(ind1) = max(min(1., icefrac(ind1)), 0.)
	1700	ENDIF
	1701	ENDIF
[3999]	1702
	1703
[5144]	1704	! ELSE ! mpc_bl_points>0
[3999]	1705
[5144]	1706	ELSE ! gaussian for environment only
[3999]	1707
[5144]	1708	alenv = (0.622 * RLVTT * qslenv(ind1)) / (rd * thetal_env(ind1)**2)
	1709	aenv = 1. / (1. + (alenv * RLVTT / rcpd))
	1710	senv = aenv * (qt_env(ind1) - qslenv(ind1))
	1711	sth = 0.
	1712	sigma1s = ratqs(ind1) * qt_env(ind1)
	1713	sigma2s = 0.
[3999]	1714
[5144]	1715	xenv = senv / (sqrt(2.) * sigma1s)
	1716	cenv(ind1) = 0.5 * (1. - 1. * erf(xenv))
	1717	ctot(ind1) = 0.5 * (1. + 1. * erf(xenv))
	1718	ctot_vol(ind1) = ctot(ind1)
	1719	qctot(ind1) = sigma1s * ((exp(-1. * xenv*2) / sqrt2pi) + xenv sqrt(2.) * cenv(ind1))
[3999]	1720
[5082]	1721	IF (ctot(ind1)<min_neb_th) THEN
[5144]	1722	ctot(ind1) = 0.
	1723	qcloud(ind1) = qslenv(ind1)
	1724	qincloud(ind1) = 0.
	1725	ELSE
	1726	qcloud(ind1) = qctot(ind1) / ctot(ind1) + qslenv(ind1)
	1727	qincloud(ind1) = MAX(qctot(ind1) / ctot(ind1), 0.)
[3999]	1728	ENDIF
	1729
[5144]	1730	ENDIF ! From the separation (thermal/envrionnement) and (environnement only,) l.335 et l.492
[3999]	1731
[5144]	1732	! Outputs used to check the PDFs
	1733	cloudth_senv(ind1) = senv
	1734	cloudth_sth(ind1) = sth
	1735	cloudth_sigmaenv(ind1) = sigma1s
	1736	cloudth_sigmath(ind1) = sigma2s
[3999]	1737
[4114]	1738	ENDDO !loop on klon
[3999]	1739
[5144]	1740	END SUBROUTINE cloudth_mpc
[4114]	1741
[3999]	1742
[5144]	1743	! ELSE ! mpc_bl_points>0
[3999]	1744
[5144]	1745	! ! Treat boundary layer mixed phase clouds
[5105]	1746
[5144]	1747	! ! thermals
	1748	! !=========
[3999]	1749
[5144]	1750	! ! ice phase
	1751	! !...........
[4910]	1752
[5144]	1753	! qiceth=0.
	1754	! deltazlev_mpc=dz(ind1,:)
	1755	! temp_mpc=ztla(ind1,:)*zpspsk(ind1,:)
	1756	! pres_mpc=pplay(ind1,:)
	1757	! fraca_mpc=fraca(ind1,:)
	1758	! snowf_mpc=snowflux(ind1,:)
	1759	! iflag_topthermals=0
	1760	! IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 0)) THEN
	1761	! iflag_topthermals = 1
	1762	! ELSE IF ((mpc_bl_points(ind1,ind2) .EQ. 1) .AND. (mpc_bl_points(ind1,ind2+1) .EQ. 1) &
	1763	! .AND. (mpc_bl_points(ind1,ind2+2) .EQ. 0) ) THEN
	1764	! iflag_topthermals = 2
	1765	! ELSE
	1766	! iflag_topthermals = 0
	1767	! ENDIF
[4910]	1768
[5144]	1769	! CALL ICE_MPC_BL_CLOUDS(ind1,ind2,klev,Ni,Ei,C_cap,d_top,iflag_topthermals,temp_mpc,pres_mpc,zqta(ind1,:), &
	1770	! qsith(ind1,:),qlth(ind1,:),deltazlev_mpc,wiceth(ind1,:),fraca_mpc,qith(ind1,:))
[5099]	1771
[5144]	1772	! ! qmax calculation
	1773	! sigma2s=(sigma2s_factor((MAX((sthl-senvl),0.)2)0.5)/((fraca(ind1,ind2)+0.02)sigma2s_power))+0.002zqta(ind1,ind2)
	1774	! deltasth=athlvert_alpha_thsigma2s
	1775	! xth1=-(sthl+deltasth)/(sqrt(2.)*sigma2s)
	1776	! xth2=-(sthl-deltasth)/(sqrt(2.)*sigma2s)
	1777	! exp_xth1 = exp(-1.xth1*2)
	1778	! exp_xth2 = exp(-1.xth2*2)
	1779	! IntJ=0.5sthl(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
	1780	! IntJ_CF=0.5(1.-1.erf(xth2))
	1781	! IntI1=(((sthl+deltasth)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
	1782	! IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp_xth1-xth2exp_xth2)
	1783	! IntI3=((sqrt2sigma2s(sthl+deltasth))/(4sqrtpideltasth))*(exp_xth1-exp_xth2)
	1784	! IntI1_CF=((sthl+deltasth)(erf(xth2)-erf(xth1)))/(4deltasth)
	1785	! IntI3_CF=(sqrt2sigma2s(exp_xth1-exp_xth2))/(4sqrtpideltasth)
	1786	! qmax=MAX(IntJ+IntI1+IntI2+IntI3,0.)
[5099]	1787
	1788
[5144]	1789	! ! Liquid phase
	1790	! !................
	1791	! ! We account for the effect of ice crystals in thermals on sthl
	1792	! ! and on the width of the distribution
[5099]	1793
[5144]	1794	! sthlc=sthl1./(1.+C_mpcqith(ind1,ind2)) &
	1795	! + (1.-1./(1.+C_mpcqith(ind1,ind2))) athl*(qsith(ind1,ind2)-qslth(ind1))
[5099]	1796
[5144]	1797	! sigma2sc=(sigma2s_factor((MAX((sthlc-senvl),0.)2)*0.5) &
	1798	! /((fraca(ind1,ind2)+0.02)**sigma2s_power)) &
	1799	! +0.002*zqta(ind1,ind2)
	1800	! deltasthc=athlvert_alpha_thsigma2sc
[5099]	1801
	1802
[5144]	1803	! xth1=-(sthlc+deltasthc)/(sqrt(2.)*sigma2sc)
	1804	! xth2=-(sthlc-deltasthc)/(sqrt(2.)*sigma2sc)
	1805	! exp_xth1 = exp(-1.xth1*2)
	1806	! exp_xth2 = exp(-1.xth2*2)
	1807	! IntJ=0.5sthlc(1-erf(xth2))+(sigma2sc/sqrt2pi)*exp_xth2
	1808	! IntJ_CF=0.5(1.-1.erf(xth2))
	1809	! IntI1=(((sthlc+deltasthc)2+(sigma2sc)2)/(8deltasthc))(erf(xth2)-erf(xth1))
	1810	! IntI2=(sigma2sc*2/(4deltasthcsqrtpi))(xth1exp_xth1-xth2exp_xth2)
	1811	! IntI3=((sqrt2sigma2sc(sthlc+deltasthc))/(4sqrtpideltasthc))*(exp_xth1-exp_xth2)
	1812	! IntI1_CF=((sthlc+deltasthc)(erf(xth2)-erf(xth1)))/(4deltasthc)
	1813	! IntI3_CF=(sqrt2sigma2sc(exp_xth1-exp_xth2))/(4sqrtpideltasthc)
	1814	! qliqth=IntJ+IntI1+IntI2+IntI3
[5099]	1815
[5144]	1816	! qlth(ind1,ind2)=MAX(0.,qliqth)
[5099]	1817
[5144]	1818	! ! Condensed water
[5099]	1819
[5144]	1820	! qcth(ind1,ind2)=qlth(ind1,ind2)+qith(ind1,ind2)
[5099]	1821
	1822
[5144]	1823	! ! consistency with subgrid distribution
[5099]	1824
[5144]	1825	! IF ((qcth(ind1,ind2) .GT. qmax) .AND. (qcth(ind1,ind2) .GT. 0)) THEN
	1826	! fraci=qith(ind1,ind2)/qcth(ind1,ind2)
	1827	! qcth(ind1,ind2)=qmax
	1828	! qith(ind1,ind2)=fraci*qmax
	1829	! qlth(ind1,ind2)=(1.-fraci)*qmax
	1830	! ENDIF
[5099]	1831
[5144]	1832	! ! Cloud Fraction
	1833	! !...............
	1834	! ! calculation of qbase which is the value of the water vapor within mixed phase clouds
	1835	! ! such that the total water in cloud = qbase+qliqth+qiceth
	1836	! ! sbase is the value of s such that int_sbase^\intfy s ds = cloud fraction
	1837	! ! sbase and qbase calculation (note that sbase is wrt liq so negative)
	1838	! ! look for an approximate solution with iteration
[5099]	1839
[5144]	1840	! ttarget=qcth(ind1,ind2)
	1841	! mini= athl*(qsith(ind1,ind2)-qslth(ind1))
	1842	! maxi= 0. !athl(3.sqrt(sigma2s))
	1843	! niter=20
	1844	! pas=(maxi-mini)/niter
	1845	! stmp=mini
	1846	! sbase=stmp
	1847	! coutref=1.E6
	1848	! DO iter=1,niter
	1849	! cout=ABS(sigma2s/SQRT(2.RPI)EXP(-((sthl-stmp)/sigma2s)*2)+(sthl-stmp)/SQRT(2.)(1.-erf(-(sthl-stmp)/sigma2s)) &
	1850	! + stmp/2.*(1.-erf(-(sthl-stmp)/sigma2s)) -ttarget)
	1851	! IF (cout .LT. coutref) THEN
	1852	! sbase=stmp
	1853	! coutref=cout
	1854	! ELSE
	1855	! stmp=stmp+pas
	1856	! ENDIF
	1857	! ENDDO
	1858	! qbase=MAX(0., sbase/athl+qslth(ind1))
[5099]	1859
[5144]	1860	! ! surface cloud fraction in thermals
	1861	! cth(ind1,ind2)=0.5*(1.-erf((sbase-sthl)/sqrt(2.)/sigma2s))
	1862	! cth(ind1,ind2)=MIN(MAX(cth(ind1,ind2),0.),1.)
[5099]	1863
	1864
[5144]	1865	! !volume cloud fraction in thermals
	1866	! !to be checked
	1867	! xth1=-(sthl+deltasth-sbase)/(sqrt(2.)*sigma2s)
	1868	! xth2=-(sthl-deltasth-sbase)/(sqrt(2.)*sigma2s)
	1869	! exp_xth1 = exp(-1.xth1*2)
	1870	! exp_xth2 = exp(-1.xth2*2)
[5099]	1871
[5144]	1872	! IntJ=0.5sthl(1-erf(xth2))+(sigma2s/sqrt2pi)*exp_xth2
	1873	! IntJ_CF=0.5(1.-1.erf(xth2))
[5099]	1874
[5144]	1875	! IF (deltasth .LT. 1.e-10) THEN
	1876	! cth_vol(ind1,ind2)=IntJ_CF
	1877	! ELSE
	1878	! IntI1=(((sthl+deltasth-sbase)2+(sigma2s)2)/(8deltasth))(erf(xth2)-erf(xth1))
	1879	! IntI2=(sigma2s*2/(4deltasthsqrtpi))(xth1exp_xth1-xth2exp_xth2)
	1880	! IntI3=((sqrt2sigma2s(sthl+deltasth))/(4sqrtpideltasth))*(exp_xth1-exp_xth2)
	1881	! IntI1_CF=((sthl-sbase+deltasth)(erf(xth2)-erf(xth1)))/(4deltasth)
	1882	! IntI3_CF=(sqrt2sigma2s(exp_xth1-exp_xth2))/(4sqrtpideltasth)
	1883	! cth_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
	1884	! ENDIF
	1885	! cth_vol(ind1,ind2)=MIN(MAX(0.,cth_vol(ind1,ind2)),1.)
[5099]	1886
	1887
	1888
[5144]	1889	! ! Environment
	1890	! !=============
	1891	! ! In the environment/downdrafts, ONLY liquid clouds
	1892	! ! See Shupe et al. 2008, JAS
[5099]	1893
[5144]	1894	! ! standard deviation of the distribution in the environment
	1895	! sth=sthl
	1896	! senv=senvl
	1897	! sigma1s_fraca = (sigma1s_factor*0.5)(fraca(ind1,ind2)**sigma1s_power) / &
	1898	! & (1-fraca(ind1,ind2))(MAX((sth-senv),0.)2)*0.5
	1899	! ! for mixed phase clouds, there is no contribution from large scale ratqs to the distribution
	1900	! ! in the environement
[5099]	1901
[5144]	1902	! sigma1s_ratqs=1E-10
	1903	! IF (cloudth_ratqsmin>0.) THEN
	1904	! sigma1s_ratqs = cloudth_ratqsmin*po(ind1)
	1905	! ENDIF
[5099]	1906
[5144]	1907	! sigma1s = sigma1s_fraca + sigma1s_ratqs
	1908	! IF (iflag_ratqs.EQ.11) THEN
	1909	! sigma1s = ratqs(ind1,ind2)po(ind1)aenv
	1910	! ENDIF
	1911	! IF (iflag_ratqs.EQ.11) THEN
	1912	! sigma1s = ratqs(ind1,ind2)po(ind1)aenvl
	1913	! ENDIF
	1914	! deltasenv=aenvlvert_alphasigma1s
	1915	! xenv1=-(senvl+deltasenv)/(sqrt(2.)*sigma1s)
	1916	! xenv2=-(senvl-deltasenv)/(sqrt(2.)*sigma1s)
	1917	! exp_xenv1 = exp(-1.xenv1*2)
	1918	! exp_xenv2 = exp(-1.xenv2*2)
[5099]	1919
[5144]	1920	! !surface CF
	1921	! cenv(ind1,ind2)=0.5(1.-1.erf(xenv1))
[5099]	1922
[5144]	1923	! !volume CF and condensed water
	1924	! IntJ=0.5senvl(1-erf(xenv2))+(sigma1s/sqrt2pi)*exp_xenv2
	1925	! IntJ_CF=0.5(1.-1.erf(xenv2))
[5099]	1926
[5144]	1927	! IF (deltasenv .LT. 1.e-10) THEN
	1928	! qcenv(ind1,ind2)=IntJ
	1929	! cenv_vol(ind1,ind2)=IntJ_CF
	1930	! ELSE
	1931	! IntI1=(((senvl+deltasenv)2+(sigma1s)2)/(8deltasenv))(erf(xenv2)-erf(xenv1))
	1932	! IntI2=(sigma1s*2/(4deltasenvsqrtpi))(xenv1exp_xenv1-xenv2exp_xenv2)
	1933	! IntI3=((sqrt2sigma1s(senv+deltasenv))/(4sqrtpideltasenv))*(exp_xenv1-exp_xenv2)
	1934	! IntI1_CF=((senvl+deltasenv)(erf(xenv2)-erf(xenv1)))/(4deltasenv)
	1935	! IntI3_CF=(sqrt2sigma1s(exp_xenv1-exp_xenv2))/(4sqrtpideltasenv)
	1936	! qcenv(ind1,ind2)=IntJ+IntI1+IntI2+IntI3 ! only liquid water in environment
	1937	! cenv_vol(ind1,ind2)=IntJ_CF+IntI1_CF+IntI3_CF
	1938	! ENDIF
[5099]	1939
[5144]	1940	! qcenv(ind1,ind2)=MAX(qcenv(ind1,ind2),0.)
	1941	! cenv_vol(ind1,ind2)=MIN(MAX(cenv_vol(ind1,ind2),0.),1.)
[5099]	1942
	1943
	1944
[5144]	1945	! ! Thermals + environment
	1946	! ! =======================
	1947	! ctot(ind1,ind2)=fraca(ind1,ind2)cth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv(ind1,ind2)
	1948	! qctot(ind1,ind2)=fraca(ind1,ind2)qcth(ind1,ind2)+(1.-1.fraca(ind1,ind2))*qcenv(ind1,ind2)
	1949	! ctot_vol(ind1,ind2)=fraca(ind1,ind2)cth_vol(ind1,ind2)+(1.-1.fraca(ind1,ind2))*cenv_vol(ind1,ind2)
	1950	! IF (qcth(ind1,ind2) .GT. 0) THEN
	1951	! icefrac(ind1,ind2)=fraca(ind1,ind2)*qith(ind1,ind2) &
	1952	! /(fraca(ind1,ind2)*qcth(ind1,ind2) &
	1953	! +(1.-1.fraca(ind1,ind2))qcenv(ind1,ind2))
	1954	! icefrac(ind1,ind2)=MAX(MIN(1.,icefrac(ind1,ind2)),0.)
	1955	! ELSE
	1956	! icefrac(ind1,ind2)=0.
	1957	! ENDIF
[5099]	1958
[5144]	1959	! IF (cenv(ind1,ind2).LT.1.e-10.OR.cth(ind1,ind2).LT.1.e-10) THEN
	1960	! ctot(ind1,ind2)=0.
	1961	! ctot_vol(ind1,ind2)=0.
	1962	! qincloud(ind1)=0.
	1963	! qcloud(ind1)=zqsatenv(ind1)
	1964	! ELSE
	1965	! qcloud(ind1)=fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)+qbase) &
	1966	! +(1.-1.fraca(ind1,ind2))(qcenv(ind1,ind2)/cenv(ind1,ind2)+qslenv(ind1))
	1967	! qincloud(ind1)=MAX(fraca(ind1,ind2)*(qcth(ind1,ind2)/cth(ind1,ind2)) &
	1968	! +(1.-1.fraca(ind1,ind2))(qcenv(ind1,ind2)/cenv(ind1,ind2)),0.)
	1969	! ENDIF
[5099]	1970
[5144]	1971	! ENDIF ! mpc_bl_points
[5099]	1972
	1973
	1974
[5144]	1975	!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
	1976	SUBROUTINE ICE_MPC_BL_CLOUDS(ind1, ind2, klev, Ni, Ei, C_cap, d_top, iflag_topthermals, temp, pres, qth, qsith, qlth, deltazlev, vith, fraca, qith)
[5099]	1977
[5144]	1978	! parameterization of ice for boundary
	1979	! layer mixed-phase clouds assuming a stationary system
[5099]	1980
[5144]	1981	! Note that vapor deposition on ice crystals and riming of liquid droplets
	1982	! depend on the ice number concentration Ni
	1983	! One could assume that Ni depends on qi, e.g., Ni=beta(qirho)**xi
	1984	! and use values from Hong et al. 2004, MWR for instance
	1985	! One may also estimate Ni as a function of T, as in Meyers 1922 or Fletcher 1962
	1986	! One could also think of a more complex expression of Ni;
	1987	! function of qi, T, the concentration in aerosols or INP ..
	1988	! Here we prefer fixing Ni to a tuning parameter
	1989	! By default we take 2.0L-1=2.0e3m-3, median value from measured vertical profiles near Svalbard
	1990	! in Mioche et al. 2017
[4910]	1991
	1992
[5144]	1993	! References:
	1994	!------------
	1995	! This parameterization is thoroughly described in Vignon et al.
[4910]	1996
[5144]	1997	! More specifically
	1998	! for the Water vapor deposition process:
[3999]	1999
[5144]	2000	! Rotstayn, L. D., 1997: A physically based scheme for the treat-
	2001	! ment of stratiform cloudfs and precipitation in large-scale
	2002	! models. I: Description and evaluation of the microphysical
	2003	! processes. Quart. J. Roy. Meteor. Soc., 123, 1227–1282.
[5099]	2004
[5144]	2005	! Morrison, H., and A. Gettelman, 2008: A new two-moment bulk
	2006	! stratiform cloud microphysics scheme in the NCAR Com-
	2007	! munity Atmosphere Model (CAM3). Part I: Description and
	2008	! numerical tests. J. Climate, 21, 3642–3659
[5099]	2009
[5144]	2010	! for the Riming process:
[5099]	2011
[5144]	2012	! Rutledge, S. A., and P. V. Hobbs, 1983: The mesoscale and micro-
	2013	! scale structure and organization of clouds and precipitation in
	2014	! midlatitude cyclones. VII: A model for the ‘‘seeder-feeder’’
	2015	! process in warm-frontal rainbands. J. Atmos. Sci., 40, 1185–1206
[5099]	2016
[5144]	2017	! Thompson, G., R. M. Rasmussen, and K. Manning, 004: Explicit
	2018	! forecasts of winter precipitation using an improved bulk
	2019	! microphysics scheme. Part I: Description and sensitivityThompson, G., R. M. Rasmussen, and K. Manning, 2004: Explicit
	2020	! forecasts of winter precipitation using an improved bulk
	2021	! microphysics scheme. Part I: Description and sensitivity analysis. Mon. Wea. Rev., 132, 519–542
[5099]	2022
[5144]	2023	! For the formation of clouds by thermals:
[5099]	2024
[5144]	2025	! Rio, C., & Hourdin, F. (2008). A thermal plume model for the convective boundary layer : Representation of cumulus clouds. Journal of
	2026	! the Atmospheric Sciences, 65, 407–425.
[5099]	2027
[5144]	2028	! Jam, A., Hourdin, F., Rio, C., & Couvreux, F. (2013). Resolved versus parametrized boundary-layer plumes. Part III: Derivation of a
	2029	! statistical scheme for cumulus clouds. Boundary-layer Meteorology, 147, 421–441. https://doi.org/10.1007/s10546-012-9789-3
[5099]	2030
	2031
	2032
[5144]	2033	! Contact: Etienne Vignon, etienne.vignon@lmd.ipsl.fr
	2034	!=============================================================================
[5099]	2035
[5117]	2036	USE phys_state_var_mod, ONLY: fm_therm, detr_therm, entr_therm
[5144]	2037	USE lmdz_yomcst
[3999]	2038
[5134]	2039	IMPLICIT NONE
[3999]	2040
[5144]	2041	INTEGER, INTENT(IN) :: ind1, ind2, klev ! horizontal and vertical indices and dimensions
[4072]	2042	INTEGER, INTENT(IN) :: iflag_topthermals ! uppermost layer of thermals ?
[5144]	2043	REAL, INTENT(IN) :: Ni ! ice number concentration [m-3]
	2044	REAL, INTENT(IN) :: Ei ! ice-droplet collision efficiency
	2045	REAL, INTENT(IN) :: C_cap ! ice crystal capacitance
	2046	REAL, INTENT(IN) :: d_top ! cloud-top ice crystal mixing parameter
	2047	REAL, DIMENSION(klev), INTENT(IN) :: temp ! temperature [K] within thermals
	2048	REAL, DIMENSION(klev), INTENT(IN) :: pres ! pressure [Pa]
	2049	REAL, DIMENSION(klev), INTENT(IN) :: qth ! mean specific water content in thermals [kg/kg]
	2050	REAL, DIMENSION(klev), INTENT(IN) :: qsith ! saturation specific humidity wrt ice in thermals [kg/kg]
	2051	REAL, DIMENSION(klev), INTENT(IN) :: qlth ! condensed liquid water in thermals, approximated value [kg/kg]
	2052	REAL, DIMENSION(klev), INTENT(IN) :: deltazlev ! layer thickness [m]
	2053	REAL, DIMENSION(klev), INTENT(IN) :: vith ! ice crystal fall velocity [m/s]
	2054	REAL, DIMENSION(klev + 1), INTENT(IN) :: fraca ! fraction of the mesh covered by thermals
	2055	REAL, DIMENSION(klev), INTENT(INOUT) :: qith ! condensed ice water , thermals [kg/kg]
[3999]	2056
[5144]	2057	INTEGER ind2p1, ind2p2
[3999]	2058	REAL rho(klev)
	2059	REAL unsurtaudet, unsurtaustardep, unsurtaurim
[4072]	2060	REAL qi, AA, BB, Ka, Dv, rhoi
	2061	REAL p0, t0, fp1, fp2
[3999]	2062	REAL alpha, flux_term
	2063	REAL det_term, precip_term, rim_term, dep_term
	2064
[5144]	2065	ind2p1 = ind2 + 1
	2066	ind2p2 = ind2 + 2
[3999]	2067
[5144]	2068	rho = pres / temp / RD ! air density kg/m3
[3999]	2069
[5144]	2070	Ka = 2.4e-2 ! thermal conductivity of the air, SI
	2071	p0 = 101325.0 ! ref pressure
	2072	T0 = 273.15 ! ref temp
	2073	rhoi = 500.0 ! cloud ice density following Reisner et al. 1998
	2074	alpha = 700. ! fallvelocity param
[3999]	2075
[5144]	2076	IF (iflag_topthermals > 0) THEN ! uppermost thermals levels
[3999]	2077
[5144]	2078	Dv = 0.0001 * 0.211 * (p0 / pres(ind2)) * ((temp(ind2) / T0)**1.94) ! water vapor diffusivity in air, SI
[3999]	2079
[5144]	2080	! Detrainment term:
[3999]	2081
[5144]	2082	unsurtaudet = detr_therm(ind1, ind2) / rho(ind2) / deltazlev(ind2)
[3999]	2083
[4072]	2084
[5144]	2085	! Deposition term
	2086	AA = RLSTT / Ka / temp(ind2) * (RLSTT / RV / temp(ind2) - 1.)
	2087	BB = 1. / (rho(ind2) * Dv * qsith(ind2))
	2088	unsurtaustardep = C_cap * (Ni*0.66) (qth(ind2) - qsith(ind2)) / qsith(ind2) &
	2089	* 4. * RPI / (AA + BB) * (6. * rho(ind2) / rhoi / RPI / Gamma(4.))**(0.33)
[3999]	2090
[5144]	2091	! Riming term neglected at this level
	2092	!unsurtaurim=rho(ind2)alpha3./rhoi/2.Eiqlth(ind2)((p0/pres(ind2))*0.4)
[3999]	2093
[5144]	2094	qi = fraca(ind2) * unsurtaustardep / MAX((d_top * unsurtaudet), 1E-12)
	2095	qi = MAX(qi, 0.)**(3. / 2.)
[3999]	2096
	2097	ELSE ! other levels, estimate qi(k) from variables at k+1 and k+2
	2098
[5144]	2099	Dv = 0.0001 * 0.211 * (p0 / pres(ind2p1)) * ((temp(ind2p1) / T0)**1.94) ! water vapor diffusivity in air, SI
[3999]	2100
[5144]	2101	! Detrainment term:
[3999]	2102
[5144]	2103	unsurtaudet = detr_therm(ind1, ind2p1) / rho(ind2p1) / deltazlev(ind2p1)
	2104	det_term = -unsurtaudet * qith(ind2p1) * rho(ind2p1)
[3999]	2105
	2106
[5144]	2107	! Deposition term
[3999]	2108
[5144]	2109	AA = RLSTT / Ka / temp(ind2p1) * (RLSTT / RV / temp(ind2p1) - 1.)
	2110	BB = 1. / (rho(ind2p1) * Dv * qsith(ind2p1))
	2111	unsurtaustardep = C_cap * (Ni*0.66) (qth(ind2p1) - qsith(ind2p1)) &
	2112	/ qsith(ind2p1) * 4. * RPI / (AA + BB) &
	2113	* (6. * rho(ind2p1) / rhoi / RPI / Gamma(4.))**(0.33)
	2114	dep_term = rho(ind2p1) * fraca(ind2p1) * (qith(ind2p1)*0.33) unsurtaustardep
[3999]	2115
[5144]	2116	! Riming term
[3999]	2117
[5144]	2118	unsurtaurim = rho(ind2p1) * alpha * 3. / rhoi / 2. * Ei * qlth(ind2p1) * ((p0 / pres(ind2p1))**0.4)
	2119	rim_term = rho(ind2p1) * fraca(ind2p1) * qith(ind2p1) * unsurtaurim
[4072]	2120
[5144]	2121	! Precip term
[3999]	2122
[5144]	2123	! We assume that there is no solid precipitation outside thermals
	2124	! so the precipitation flux within thermals is equal to the precipitation flux
	2125	! at mesh-scale divided by thermals fraction
[3999]	2126
[5144]	2127	fp2 = 0.
	2128	fp1 = 0.
	2129	IF (fraca(ind2p1) > 0.) THEN
	2130	fp2 = -qith(ind2p2) * rho(ind2p2) * vith(ind2p2) * fraca(ind2p2)! flux defined positive upward
	2131	fp1 = -qith(ind2p1) * rho(ind2p1) * vith(ind2p1) * fraca(ind2p1)
	2132	ENDIF
[3999]	2133
[5144]	2134	precip_term = -1. / deltazlev(ind2p1) * (fp2 - fp1)
[4072]	2135
[5144]	2136	! Calculation in a top-to-bottom loop
[3999]	2137
[5144]	2138	IF (fm_therm(ind1, ind2p1) > 0.) THEN
	2139	qi = 1. / fm_therm(ind1, ind2p1) * &
	2140	(deltazlev(ind2p1) * (-rim_term - dep_term - det_term - precip_term) + &
	2141	fm_therm(ind1, ind2p2) * (qith(ind2p1)))
	2142	END IF
[5105]	2143
[5144]	2144	ENDIF ! top thermals
[3999]	2145
[5144]	2146	qith(ind2) = MAX(0., qi)
	2147
	2148	END SUBROUTINE ICE_MPC_BL_CLOUDS
	2149
[4651]	2150	END MODULE lmdz_cloudth

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

Download in other formats:

Original Format