Context Navigation

fisrtilp.F90 @ 2072

Last change on this file since 2072 was 2006, checked in by fhourdin, 11 years ago

Modification de la spécification de la plage de température pour
la phase mixte liquide/glace des nuages.
Contrôle par les paramètres t_glace_min/max, exposant_glace, iflag_t_glace

Modifying the specification of the mixte liquid/ice phase for cloud water.

Jean-Baptiste Madeleine

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 33.3 KB

Rev	Line
[524]	1	!
[1403]	2	! $Id: fisrtilp.F90 2006 2014-04-06 11:31:28Z lguez $
[524]	3	!
[1472]	4	!
	5	SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs, &
[1742]	6	d_t, d_q, d_ql, rneb, radliq, rain, snow, &
	7	pfrac_impa, pfrac_nucl, pfrac_1nucl, &
	8	frac_impa, frac_nucl, beta, &
	9	prfl, psfl, rhcl, zqta, fraca, &
[1849]	10	ztv, zpspsk, ztla, zthl, iflag_cldcon, &
	11	iflag_ice_thermo)
[524]	12
[1472]	13	!
	14	USE dimphy
[2006]	15	USE microphys_mod ! cloud microphysics (JBM 3/14)
[1472]	16	IMPLICIT none
	17	!======================================================================
	18	! Auteur(s): Z.X. Li (LMD/CNRS)
	19	! Date: le 20 mars 1995
	20	! Objet: condensation et precipitation stratiforme.
	21	! schema de nuage
	22	!======================================================================
	23	!======================================================================
	24	!ym include "dimensions.h"
	25	!ym include "dimphy.h"
	26	include "YOMCST.h"
	27	include "tracstoke.h"
	28	include "fisrtilp.h"
[2006]	29	include "nuage.h" ! JBM (3/14)
[1506]	30	include "iniprint.h"
	31
[1472]	32	!
	33	! Arguments:
	34	!
	35	REAL dtime ! intervalle du temps (s)
	36	REAL paprs(klon,klev+1) ! pression a inter-couche
	37	REAL pplay(klon,klev) ! pression au milieu de couche
	38	REAL t(klon,klev) ! temperature (K)
	39	REAL q(klon,klev) ! humidite specifique (kg/kg)
	40	REAL d_t(klon,klev) ! incrementation de la temperature (K)
	41	REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
	42	REAL d_ql(klon,klev) ! incrementation de l'eau liquide
	43	REAL rneb(klon,klev) ! fraction nuageuse
	44	REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements
	45	REAL rhcl(klon,klev) ! humidite relative en ciel clair
	46	REAL rain(klon) ! pluies (mm/s)
	47	REAL snow(klon) ! neige (mm/s)
	48	REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
	49	REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
	50	REAL ztv(klon,klev)
	51	REAL zqta(klon,klev),fraca(klon,klev)
	52	REAL sigma1(klon,klev),sigma2(klon,klev)
	53	REAL qltot(klon,klev),ctot(klon,klev)
	54	REAL zpspsk(klon,klev),ztla(klon,klev)
	55	REAL zthl(klon,klev)
[1849]	56	REAL ztfondue, qsl, qsi
[1403]	57
[1472]	58	logical lognormale(klon)
[1849]	59	logical ice_thermo
[1411]	60
[1472]	61	!AA
	62	! Coeffients de fraction lessivee : pour OFF-LINE
	63	!
	64	REAL pfrac_nucl(klon,klev)
	65	REAL pfrac_1nucl(klon,klev)
	66	REAL pfrac_impa(klon,klev)
	67	!
	68	! Fraction d'aerosols lessivee par impaction et par nucleation
	69	! POur ON-LINE
	70	!
	71	REAL frac_impa(klon,klev)
	72	REAL frac_nucl(klon,klev)
	73	real zct ,zcl
	74	!AA
	75	!
	76	! Options du programme:
	77	!
	78	REAL seuil_neb ! un nuage existe vraiment au-dela
	79	PARAMETER (seuil_neb=0.001)
[524]	80
[1472]	81	INTEGER ninter ! sous-intervals pour la precipitation
	82	INTEGER ncoreczq
	83	INTEGER iflag_cldcon
[1849]	84	INTEGER iflag_ice_thermo
[1472]	85	PARAMETER (ninter=5)
	86	LOGICAL evap_prec ! evaporation de la pluie
	87	PARAMETER (evap_prec=.TRUE.)
	88	REAL ratqs(klon,klev) ! determine la largeur de distribution de vapeur
	89	logical ptconv(klon,klev) ! determine la largeur de distribution de vapeur
[524]	90
[1472]	91	real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon)
	92	real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon)
	93	real erf
	94	REAL qcloud(klon)
	95	!
	96	LOGICAL cpartiel ! condensation partielle
	97	PARAMETER (cpartiel=.TRUE.)
	98	REAL t_coup
	99	PARAMETER (t_coup=234.0)
	100	!
	101	! Variables locales:
	102	!
	103	INTEGER i, k, n, kk
[1901]	104	REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
	105	REAL Tbef(klon),qlbef(klon),DT(klon),num,denom
	106	LOGICAL convergence(klon)
	107	REAL DDT0
	108	PARAMETER (DDT0=.01)
	109	INTEGER n_i, iter
	110
[1849]	111	REAL zrfl(klon), zrfln(klon), zqev, zqevt
	112	REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti
	113	REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
	114	REAL zoliqp(klon), zoliqi(klon)
[2006]	115	REAL zt(klon)
	116	! JBM (3/14) nexpo is replaced by exposant_glace
	117	! REAL nexpo ! exponentiel pour glace/eau
	118	! INTEGER, PARAMETER :: nexpo=6
	119	INTEGER exposant_glace_old
	120	REAL t_glace_min_old
[1472]	121	REAL zdz(klon),zrho(klon),ztot , zrhol(klon)
	122	REAL zchau ,zfroi ,zfice(klon),zneb(klon)
[1849]	123	REAL zmelt, zpluie, zice, zcondold
	124	PARAMETER (ztfondue=278.15)
[1472]	125	!
	126	LOGICAL appel1er
	127	SAVE appel1er
	128	!$OMP THREADPRIVATE(appel1er)
	129	!
	130	!---------------------------------------------------------------
	131	!
	132	!AA Variables traceurs:
	133	!AA Provisoire !!! Parametres alpha du lessivage
	134	!AA A priori on a 4 scavenging # possibles
	135	!
	136	REAL a_tr_sca(4)
	137	save a_tr_sca
	138	!$OMP THREADPRIVATE(a_tr_sca)
	139	!
	140	! Variables intermediaires
	141	!
	142	REAL zalpha_tr
	143	REAL zfrac_lessi
	144	REAL zprec_cond(klon)
	145	!AA
[1742]	146	! RomP >>> 15 nov 2012
	147	REAL beta(klon,klev) ! taux de conversion de l'eau cond
	148	! RomP <<<
[1472]	149	REAL zmair, zcpair, zcpeau
	150	! Pour la conversion eau-neige
	151	REAL zlh_solid(klon), zm_solid
	152	!IM
	153	!ym INTEGER klevm1
	154	!---------------------------------------------------------------
	155	!
	156	! Fonctions en ligne:
	157	!
	158	REAL fallvs,fallvc ! vitesse de chute pour crystaux de glace
	159	REAL zzz
	160	include "YOETHF.h"
	161	include "FCTTRE.h"
	162	fallvc (zzz) = 3.29/2.0 * ((zzz)*0.16) ffallv_con
	163	fallvs (zzz) = 3.29/2.0 * ((zzz)*0.16) ffallv_lsc
	164	!
	165	DATA appel1er /.TRUE./
	166	!ym
[1849]	167	ice_thermo = iflag_ice_thermo .GE. 1
[1472]	168	zdelq=0.0
[524]	169
[1506]	170	if (prt_level>9)write(lunout,*)'NUAGES4 A. JAM'
[1472]	171	IF (appel1er) THEN
	172	!
[1575]	173	WRITE(lunout,*) 'fisrtilp, ninter:', ninter
	174	WRITE(lunout,*) 'fisrtilp, evap_prec:', evap_prec
	175	WRITE(lunout,*) 'fisrtilp, cpartiel:', cpartiel
[1472]	176	IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN
[1575]	177	WRITE(lunout,*) 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
	178	WRITE(lunout,*) 'Je prefere un sous-intervalle de 6 minutes'
[1472]	179	! CALL abort
	180	ENDIF
	181	appel1er = .FALSE.
	182	!
	183	!AA initialiation provisoire
	184	a_tr_sca(1) = -0.5
	185	a_tr_sca(2) = -0.5
	186	a_tr_sca(3) = -0.5
	187	a_tr_sca(4) = -0.5
	188	!
	189	!AA Initialisation a 1 des coefs des fractions lessivees
	190	!
	191	!cdir collapse
	192	DO k = 1, klev
	193	DO i = 1, klon
	194	pfrac_nucl(i,k)=1.
	195	pfrac_1nucl(i,k)=1.
	196	pfrac_impa(i,k)=1.
[1742]	197	beta(i,k)=0. !RomP initialisation
[1472]	198	ENDDO
	199	ENDDO
[524]	200
[1472]	201	ENDIF ! test sur appel1er
	202	!
	203	!MAf Initialisation a 0 de zoliq
	204	! DO i = 1, klon
	205	! zoliq(i)=0.
	206	! ENDDO
	207	! Determiner les nuages froids par leur temperature
	208	! nexpo regle la raideur de la transition eau liquide / eau glace.
	209	!
[2006]	210	IF (iflag_t_glace.EQ.0) THEN
	211	! ztglace = RTT - 15.0
	212	t_glace_min_old = RTT - 15.0
	213	!AJ<
	214	IF (ice_thermo) THEN
	215	! nexpo = 2
	216	exposant_glace_old = 2
	217	ELSE
	218	! nexpo = 6
	219	exposant_glace_old = 6
	220	ENDIF
[1849]	221	ENDIF
[2006]	222
[1849]	223	!! RLVTT = 2.501e6 ! pas de redefinition des constantes physiques (jyg)
	224	!! RLSTT = 2.834e6 ! pas de redefinition des constantes physiques (jyg)
	225	!>AJ
[1472]	226	!cc nexpo = 1
	227	!
	228	! Initialiser les sorties:
	229	!
	230	!cdir collapse
	231	DO k = 1, klev+1
	232	DO i = 1, klon
	233	prfl(i,k) = 0.0
	234	psfl(i,k) = 0.0
	235	ENDDO
	236	ENDDO
[524]	237
[1472]	238	!cdir collapse
	239	DO k = 1, klev
	240	DO i = 1, klon
	241	d_t(i,k) = 0.0
	242	d_q(i,k) = 0.0
	243	d_ql(i,k) = 0.0
	244	rneb(i,k) = 0.0
	245	radliq(i,k) = 0.0
	246	frac_nucl(i,k) = 1.
	247	frac_impa(i,k) = 1.
	248	ENDDO
	249	ENDDO
	250	DO i = 1, klon
	251	rain(i) = 0.0
	252	snow(i) = 0.0
	253	zoliq(i)=0.
	254	! ENDDO
	255	!
	256	! Initialiser le flux de precipitation a zero
	257	!
	258	! DO i = 1, klon
	259	zrfl(i) = 0.0
[1849]	260	zifl(i) = 0.0
[1472]	261	zneb(i) = seuil_neb
	262	ENDDO
	263	!
	264	!
	265	!AA Pour plus de securite
[524]	266
[1472]	267	zalpha_tr = 0.
	268	zfrac_lessi = 0.
[524]	269
[1472]	270	!AA----------------------------------------------------------
	271	!
	272	ncoreczq=0
	273	! Boucle verticale (du haut vers le bas)
	274	!
	275	!IM : klevm1
	276	!ym klevm1=klev-1
	277	DO k = klev, 1, -1
	278	!
	279	!AA----------------------------------------------------------
	280	!
	281	DO i = 1, klon
	282	zt(i)=t(i,k)
	283	zq(i)=q(i,k)
	284	ENDDO
	285	!
	286	! Calculer la varition de temp. de l'air du a la chaleur sensible
	287	! transporter par la pluie.
	288	! Il resterait a rajouter cet effet de la chaleur sensible sur les
	289	! flux de surface, du a la diff. de temp. entre le 1er niveau et la
	290	! surface.
	291	!
	292	IF(k.LE.klevm1) THEN
	293	DO i = 1, klon
	294	!IM
	295	zmair=(paprs(i,k)-paprs(i,k+1))/RG
	296	zcpair=RCPD(1.0+RVTMP2zq(i))
	297	zcpeau=RCPD*RVTMP2
	298	zt(i) = ( (t(i,k+1)+d_t(i,k+1))zrfl(i)dtime*zcpeau &
	299	+ zmairzcpairzt(i) ) &
	300	/ (zmairzcpair + zrfl(i)dtime*zcpeau)
	301	! C WRITE (6,*) 'cppluie ', zt(i)-(t(i,k+1)+d_t(i,k+1))
	302	ENDDO
	303	ENDIF
	304	!
	305	!
	306	! Calculer l'evaporation de la precipitation
	307	!
[524]	308
	309
[1849]	310	! Calculer l'evaporation de la precipitation
	311	!
	312
	313
[1472]	314	IF (evap_prec) THEN
	315	DO i = 1, klon
[1849]	316	!AJ<
	317	!! IF (zrfl(i) .GT.0.) THEN
	318	IF (zrfl(i)+zifl(i).GT.0.) THEN
	319	!>AJ
[1472]	320	IF (thermcep) THEN
	321	zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
	322	zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
	323	zqs(i)=MIN(0.5,zqs(i))
	324	zcor=1./(1.-RETV*zqs(i))
	325	zqs(i)=zqs(i)*zcor
	326	ELSE
	327	IF (zt(i) .LT. t_coup) THEN
	328	zqs(i) = qsats(zt(i)) / pplay(i,k)
	329	ELSE
	330	zqs(i) = qsatl(zt(i)) / pplay(i,k)
	331	ENDIF
	332	ENDIF
[1849]	333	ENDIF ! (zrfl(i)+zifl(i).GT.0.)
	334	ENDDO
	335	!AJ<
	336	IF (.NOT. ice_thermo) THEN
	337	DO i = 1, klon
	338	!AJ<
	339	!! IF (zrfl(i) .GT.0.) THEN
	340	IF (zrfl(i)+zifl(i).GT.0.) THEN
	341	!>AJ
	342	zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
	343	zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) &
	344	* (paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)RD/RG
	345	zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
	346	* RG*dtime/(paprs(i,k)-paprs(i,k+1))
	347	zqev = MIN (zqev, zqevt)
	348	zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
	349	/RG/dtime
	350
	351	! pour la glace, on ré-évapore toute la précip dans la
	352	! couche du dessous
	353	! la glace venant de la couche du dessus est simplement
	354	! dans la couche du dessous.
	355
	356	IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
	357
	358	zq(i) = zq(i) - (zrfln(i)-zrfl(i)) &
	359	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
	360	zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
	361	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
	362	* RLVTT/RCPD/(1.0+RVTMP2*zq(i))
	363	zrfl(i) = zrfln(i)
	364	zifl(i) = 0.
	365	ENDIF ! (zrfl(i)+zifl(i).GT.0.)
	366	ENDDO
	367	!
	368	ELSE ! (.NOT. ice_thermo)
	369	!
	370	DO i = 1, klon
	371	!AJ<
	372	!! IF (zrfl(i) .GT.0.) THEN
	373	IF (zrfl(i)+zifl(i).GT.0.) THEN
	374	!>AJ
	375	!JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	376	! Modification de l'évaporation avec la glace
	377	! Différentiation entre précipitation liquide et solide
	378	! On suppose que coef_evai=2*coef_eva
	379	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	380
	381	zqev0 = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
	382	! zqev0 = MAX (0.0, zqs(i)-zq(i) )
[524]	383
[1849]	384	!JAM !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	385	! On différencie qsat pour l'eau et la glace
	386	! Si zdelta=1. --> glace
	387	! Si zdelta=0. --> eau liquide
	388	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	389
	390	qsl= R2ES*FOEEW(zt(i),0.)/pplay(i,k)
	391	qsl= MIN(0.5,qsl)
	392	zcor= 1./(1.-RETV*qsl)
	393	qsl= qsl*zcor
	394
	395	zqevt = 1.coef_eva(1.0-zq(i)/qsl)*SQRT(zrfl(i)) &
	396	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
	397	zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) &
	398	RGdtime/(paprs(i,k)-paprs(i,k+1))
[524]	399
[1849]	400	qsi= R2ES*FOEEW(zt(i),1.)/pplay(i,k)
	401	qsi= MIN(0.5,qsi)
	402	zcor= 1./(1.-RETV*qsi)
	403	qsi= qsi*zcor
[1472]	404
[1849]	405	zqevti = 1.coef_eva(1.0-zq(i)/qsi)*SQRT(zifl(i)) &
	406	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
	407	zqevti = MAX(0.0,MIN(zqevti,zifl(i))) &
	408	RGdtime/(paprs(i,k)-paprs(i,k+1))
	409
	410
	411	!JAM!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	412	! Vérification sur l'évaporation
	413	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	414
	415	IF (zqevt+zqevti.GT.zqev0) THEN
	416	zqev=zqev0*zqevt/(zqevt+zqevti)
	417	zqevi=zqev0*zqevti/(zqevt+zqevti)
	418
	419	ELSE
	420	IF (zqevt+zqevti.GT.0.) THEN
	421	zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt)
	422	zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti)
	423	ELSE
	424	zqev=0.
	425	zqevi=0.
	426	ENDIF
	427	ENDIF
	428
	429	zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) &
	430	/RG/dtime)
	431	zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) &
	432	/RG/dtime)
	433
	434	! Pour la glace, on révapore toute la précip dans la couche du dessous
	435	! la glace venant de la couche du dessus est simplement dans la couche
	436	! du dessous.
	437
	438	! IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0.
	439	! print*,zrfl(i),zrfln(i),zqevt,zqevti,RLMLT,'fluxdeprecip'
	440	zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) &
	441	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
	442	zt(i) = zt(i) + (zrfln(i)-zrfl(i)) &
	443	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
	444	* RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
	445	+ (zifln(i)-zifl(i)) &
	446	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
	447	* RLSTT/RCPD/(1.0+RVTMP2*zq(i))
	448
	449	zrfl(i) = zrfln(i)
	450	zifl(i) = zifln(i)
	451
	452	ENDIF ! (zrfl(i)+zifl(i).GT.0.)
[1472]	453	ENDDO
[1849]	454
	455	ENDIF ! (.NOT. ice_thermo)
	456
	457	ENDIF ! (evap_prec)
[1472]	458	!
	459	! Calculer Qs et L/Cp*dQs/dT:
	460	!
	461	IF (thermcep) THEN
	462	DO i = 1, klon
[524]	463	zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
	464	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
	465	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
	466	zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
	467	zqs(i) = MIN(0.5,zqs(i))
	468	zcor = 1./(1.-RETV*zqs(i))
	469	zqs(i) = zqs(i)*zcor
	470	zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
[1472]	471	ENDDO
	472	ELSE
	473	DO i = 1, klon
	474	IF (zt(i).LT.t_coup) THEN
	475	zqs(i) = qsats(zt(i))/pplay(i,k)
	476	zdqs(i) = dqsats(zt(i),zqs(i))
	477	ELSE
	478	zqs(i) = qsatl(zt(i))/pplay(i,k)
	479	zdqs(i) = dqsatl(zt(i),zqs(i))
	480	ENDIF
	481	ENDDO
	482	ENDIF
	483	!
	484	! Determiner la condensation partielle et calculer la quantite
	485	! de l'eau condensee:
	486	!
[1901]	487	!verification de la valeur de iflag_fisrtilp_qsat pour iflag_ice_thermo=1
	488	if ((iflag_ice_thermo.eq.1).and.(iflag_fisrtilp_qsat.ne.0)) then
	489	write(,) " iflag_ice_thermo==1 requires iflag_fisrtilp_qsat==0", &
	490	" but iflag_fisrtilp_qsat=",iflag_fisrtilp_qsat, ". Might as well stop here."
	491	stop
	492	endif
[1403]	493
[1472]	494	IF (cpartiel) THEN
[524]	495
[1472]	496	! print*,'Dans partiel k=',k
	497	!
	498	! Calcul de l'eau condensee et de la fraction nuageuse et de l'eau
	499	! nuageuse a partir des PDF de Sandrine Bony.
	500	! rneb : fraction nuageuse
	501	! zqn : eau totale dans le nuage
	502	! zcond : eau condensee moyenne dans la maille.
	503	! on prend en compte le réchauffement qui diminue la partie
	504	! condensee
	505	!
	506	! Version avec les raqts
[524]	507
[1472]	508	if (iflag_pdf.eq.0) then
[524]	509
	510	do i=1,klon
[1472]	511	zdelq = min(ratqs(i,k),0.99) * zq(i)
	512	rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
	513	zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
[524]	514	enddo
	515
[1472]	516	else
	517	!
	518	! Version avec les nouvelles PDFs.
[524]	519	do i=1,klon
	520	if(zq(i).lt.1.e-15) then
[1472]	521	ncoreczq=ncoreczq+1
	522	zq(i)=1.e-15
[524]	523	endif
[1472]	524	enddo
[1403]	525
[1472]	526	if (iflag_cldcon>=5) then
[1403]	527
[1472]	528	call cloudth(klon,klev,k,ztv, &
	529	zq,zqta,fraca, &
	530	qcloud,ctot,zpspsk,paprs,ztla,zthl, &
	531	ratqs,zqs,t)
[1403]	532
[1472]	533	do i=1,klon
[1403]	534	rneb(i,k)=ctot(i,k)
	535	zqn(i)=qcloud(i)
[1472]	536	enddo
[1403]	537
[1472]	538	endif
	539
	540	if (iflag_cldcon <= 4) then
	541	lognormale = .true.
[1507]	542	elseif (iflag_cldcon >= 6) then
[1472]	543	! lognormale en l'absence des thermiques
	544	lognormale = fraca(:,k) < 1e-10
	545	else
	546	! Dans le cas iflag_cldcon=5, on prend systématiquement la
	547	! bi-gaussienne
	548	lognormale = .false.
	549	end if
	550
	551	do i=1,klon
[1901]	552	Tbef(i)=zt(i)
	553	qlbef(i)=0.
[1472]	554	if (lognormale(i)) then
	555	zpdf_sig(i)=ratqs(i,k)*zq(i)
	556	zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
	557	zpdf_delta(i)=log(zq(i)/zqs(i))
	558	zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
	559	zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
	560	zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
	561	zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
	562	zpdf_e1(i)=1.-erf(zpdf_e1(i))
	563	zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
	564	zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
	565	zpdf_e2(i)=1.-erf(zpdf_e2(i))
[1901]	566
	567	if (zpdf_e1(i).lt.1.e-10) then
	568	rneb(i,k)=0.
	569	zqn(i)=zqs(i)
	570	else
	571	rneb(i,k)=0.5*zpdf_e1(i)
	572	zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
	573	endif
	574
	575	qlbef(i)=max(0.,zqn(i)-zqs(i))
	576	num=-Tbef(i)+zt(i)+rneb(i,k)RLVTT/RCPD/(1.0+RVTMP2zq(i))*qlbef(i)
	577	denom=1.+rneb(i,k)*zdqs(i)
	578	DT(i)=num/denom
[1411]	579	endif
[1472]	580	enddo
[1901]	581
	582	n_i=1
	583	! do while ((abs(DT(i)).gt.DDT0).and.((zqn(i)-zqs(i)).gt.0.))
	584	! iflag_fisrtilp_qsat=0: qsat ne varie pas avec T
	585	! iflag_fisrtilp_qsat > 1 : nombre d iterations max pour converger sur le calcul de qsat(T)
	586	! do while (n_i.le.iflag_fisrtilp_qsat)
	587	if (iflag_fisrtilp_qsat.ge.1) then
	588	do iter=1,iflag_fisrtilp_qsat
	589	do i=1,klon
	590	convergence(i)=abs(DT(i)).gt.DDT0
	591	if (convergence(i).and.lognormale(i)) then
	592	Tbef(i)=Tbef(i)+DT(i)
	593	zdelta=MAX(0.,SIGN(1.,RTT-Tbef(i)))
	594	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
	595	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
	596	zqs(i)= R2ES * FOEEW(Tbef(i),zdelta)/pplay(i,k)
	597	zqs(i)=MIN(0.5,zqs(i))
	598	zcor=1./(1.-retv*zqs(i))
	599	zqs(i)=zqs(i)*zcor
	600	zdqs(i)=FOEDE(Tbef(i),zdelta,zcvm5,zqs(i),zcor)
[1403]	601
[1901]	602	zpdf_sig(i)=ratqs(i,k)*zq(i)
	603	zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2))
	604	zpdf_delta(i)=log(zq(i)/zqs(i))
	605	zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.))
	606	zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.))
	607	zpdf_e1(i)=zpdf_a(i)-zpdf_b(i)
	608	zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i))
	609	zpdf_e1(i)=1.-erf(zpdf_e1(i))
	610	zpdf_e2(i)=zpdf_a(i)+zpdf_b(i)
	611	zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i))
	612	zpdf_e2(i)=1.-erf(zpdf_e2(i))
	613
[1472]	614	if (zpdf_e1(i).lt.1.e-10) then
	615	rneb(i,k)=0.
	616	zqn(i)=zqs(i)
	617	else
	618	rneb(i,k)=0.5*zpdf_e1(i)
	619	zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i)
	620	endif
[1411]	621
[1901]	622	qlbef(i)=max(0.,zqn(i)-zqs(i))
	623	num=-Tbef(i)+zt(i)+rneb(i,k)RLVTT/RCPD/(1.0+RVTMP2zq(i))*qlbef(i)
	624	denom=1.+rneb(i,k)*zdqs(i)
	625	DT(i)=num/denom
	626	n_i=n_i+1
	627	endif
	628	enddo
	629	enddo
[1411]	630
[1901]	631	endif
[524]	632
[1901]	633
[524]	634	endif ! iflag_pdf
	635
[1901]	636	if (iflag_fisrtilp_qsat.eq.-1) then
	637	!CR: ATTENTION option fausse mais a existe
[1146]	638	DO i=1,klon
	639	IF (rneb(i,k) .LE. 0.0) THEN
	640	zqn(i) = 0.0
	641	rneb(i,k) = 0.0
	642	zcond(i) = 0.0
	643	rhcl(i,k)=zq(i)/zqs(i)
	644	ELSE IF (rneb(i,k) .GE. 1.0) THEN
	645	zqn(i) = zq(i)
[1901]	646	rneb(i,k) = 1.0
	647	zcond(i) = MAX(0.0,zqn(i)-zqs(i))/(1+zdqs(i))
[1146]	648	rhcl(i,k)=1.0
	649	ELSE
[1901]	650	zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1+zdqs(i))
[1146]	651	rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
	652	ENDIF
	653	ENDDO
[1901]	654
	655	ELSE
	656
	657	DO i=1,klon
	658	IF (rneb(i,k) .LE. 0.0) THEN
	659	zqn(i) = 0.0
	660	rneb(i,k) = 0.0
	661	zcond(i) = 0.0
	662	rhcl(i,k)=zq(i)/zqs(i)
	663	ELSE IF (rneb(i,k) .GE. 1.0) THEN
	664	zqn(i) = zq(i)
	665	rneb(i,k) = 1.0
	666	zcond(i) = MAX(0.0,zqn(i)-zqs(i))
	667	rhcl(i,k)=1.0
	668	ELSE
	669	zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
	670	rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
	671	ENDIF
	672	ENDDO
	673
	674
	675	ENDIF
	676
[1472]	677	! do i=1,klon
	678	! IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
	679	! IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
	680	! rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
	681	!c zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
	682	!c On ne divise pas par 1+zdqs pour forcer a avoir l'eau predite par
	683	!c la convection.
	684	!c ATTENTION !!! Il va falloir verifier tout ca.
	685	! zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)
	686	!c print*,'ZDQS ',zdqs(i)
	687	!c--Olivier
	688	! rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i)
	689	! IF (rneb(i,k) .LE. 0.0) rhcl(i,k)=zq(i)/zqs(i)
	690	! IF (rneb(i,k) .GE. 1.0) rhcl(i,k)=1.0
	691	!c--fin
	692	! ENDDO
	693	ELSE
	694	DO i = 1, klon
	695	IF (zq(i).GT.zqs(i)) THEN
	696	rneb(i,k) = 1.0
	697	ELSE
	698	rneb(i,k) = 0.0
	699	ENDIF
	700	zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
	701	ENDDO
	702	ENDIF
	703	!
	704	DO i = 1, klon
	705	zq(i) = zq(i) - zcond(i)
	706	! zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
	707	ENDDO
[1849]	708	!AJ<
	709	IF (.NOT. ice_thermo) THEN
[1901]	710	if (iflag_fisrtilp_qsat.lt.1) then
	711	DO i = 1, klon
	712	zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
	713	ENDDO
	714	else if (iflag_fisrtilp_qsat.gt.0) then
	715	DO i= 1, klon
	716	if (lognormale(i)) then
	717	zt(i)=Tbef(i)
	718	else
	719	zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
	720	endif
	721	ENDDO
	722	endif
[1849]	723	ELSE
[2006]	724	IF (iflag_t_glace.EQ.0) THEN
[1901]	725	if (iflag_fisrtilp_qsat.lt.1) then
	726	DO i = 1, klon
[2006]	727	zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (273.15-t_glace_min_old)
[1901]	728	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
[2006]	729	zfice(i) = zfice(i)**exposant_glace_old
	730	! zfice(i) = zfice(i)**nexpo
[1901]	731	zt(i) = zt(i) + (1.-zfice(i))zcond(i) RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
[1849]	732	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*zq(i))
[1901]	733	ENDDO
	734	else
	735	DO i=1, klon
[2006]	736	zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (273.15-t_glace_min_old)
[1901]	737	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
[2006]	738	zfice(i) = zfice(i)**exposant_glace_old
	739	! zfice(i) = zfice(i)**nexpo
[1901]	740	!CR: ATTENTION zt different de Tbef: à corriger
	741	zt(i) = zt(i) + (1.-zfice(i))zcond(i) RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i))) &
	742	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
	743	ENDDO
	744	endif
[1849]	745	! print*,zt(i),zrfl(i),zifl(i),'temp1'
[2006]	746	ELSE ! of IF (iflag_t_glace.EQ.0)
	747	if (iflag_fisrtilp_qsat.lt.1) then
	748	DO i = 1, klon
	749	! JBM: icefrac_lsc is now a function contained in microphys_mod
	750	zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
	751	t_glace_max, exposant_glace)
	752	zt(i) = zt(i) + (1.-zfice(i))zcond(i) RLVTT/RCPD/(1.0+RVTMP2*zq(i)) &
	753	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*zq(i))
	754	ENDDO
	755	else
	756	DO i=1, klon
	757	! JBM: icefrac_lsc is now a function contained in microphys_mod
	758	zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
	759	t_glace_max, exposant_glace)
	760	!CR: ATTENTION zt different de Tbef: à corriger
	761	zt(i) = zt(i) + (1.-zfice(i))zcond(i) RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i))) &
	762	+zfice(i)zcond(i) RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i)))
	763	ENDDO
	764	endif
	765	! print*,zt(i),zrfl(i),zifl(i),'temp1'
	766	ENDIF
[1849]	767	ENDIF
	768	!>AJ
[1472]	769	!
	770	! Partager l'eau condensee en precipitation et eau liquide nuageuse
	771	!
	772	DO i = 1, klon
	773	IF (rneb(i,k).GT.0.0) THEN
	774	zoliq(i) = zcond(i)
	775	zrho(i) = pplay(i,k) / zt(i) / RD
	776	zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
[1849]	777	ENDIF
	778	ENDDO
	779	!AJ<
	780	IF (.NOT. ice_thermo) THEN
[2006]	781	IF (iflag_t_glace.EQ.0) THEN
	782	DO i = 1, klon
	783	IF (rneb(i,k).GT.0.0) THEN
	784	zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (273.13-t_glace_min_old)
	785	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
	786	zfice(i) = zfice(i)**exposant_glace_old
	787	! zfice(i) = zfice(i)**nexpo
	788	!! zfice(i)=0.
	789	ENDIF
	790	ENDDO
	791	ELSE ! of IF (iflag_t_glace.EQ.0)
	792	DO i = 1, klon
	793	IF (rneb(i,k).GT.0.0) THEN
	794	! JBM: icefrac_lsc is now a function contained in microphys_mod
	795	zfice(i) = icefrac_lsc(zt(i), t_glace_min, &
	796	t_glace_max, exposant_glace)
	797	ENDIF
	798	ENDDO
	799	ENDIF
[1849]	800	ENDIF
	801	DO i = 1, klon
	802	IF (rneb(i,k).GT.0.0) THEN
[1472]	803	zneb(i) = MAX(rneb(i,k), seuil_neb)
[1849]	804	! zt(i) = zt(i)+zcond(i)zfice(i)RLMLT/RCPD/(1.0+RVTMP2*zq(i))
	805	! print*,zt(i),'fractionglace'
	806	!>AJ
[1472]	807	radliq(i,k) = zoliq(i)/REAL(ninter+1)
	808	ENDIF
	809	ENDDO
	810	!
	811	DO n = 1, ninter
	812	DO i = 1, klon
	813	IF (rneb(i,k).GT.0.0) THEN
	814	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
[1855]	815	! Initialization of zpluie and zice:
	816	zpluie=0
	817	zice=0
[1472]	818	IF (zneb(i).EQ.seuil_neb) THEN
	819	ztot = 0.0
	820	ELSE
	821	! quantite d'eau a eliminer: zchau
	822	! meme chose pour la glace: zfroi
	823	if (ptconv(i,k)) then
	824	zcl =cld_lc_con
	825	zct =1./cld_tau_con
	826	zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
	827	fallvc(zrhol(i)) zfice(i)
	828	else
	829	zcl =cld_lc_lsc
	830	zct =1./cld_tau_lsc
	831	zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) &
	832	fallvs(zrhol(i)) zfice(i)
	833	endif
	834	zchau = zct dtime/REAL(ninter) zoliq(i) &
	835	(1.0-EXP(-(zoliq(i)/zneb(i)/zcl )2)) (1.-zfice(i))
[1849]	836	!AJ<
	837	IF (.NOT. ice_thermo) THEN
	838	ztot = zchau + zfroi
	839	ELSE
	840	zpluie = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i)))
	841	zice = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i))
	842	ztot = zpluie + zice
	843	ENDIF
	844	!>AJ
[1472]	845	ztot = MAX(ztot ,0.0)
	846	ENDIF
	847	ztot = MIN(ztot,zoliq(i))
[1849]	848	!AJ<
	849	! zoliqp = MAX(zoliq(i)(1.-zfice(i))-1.zpluie , 0.0)
	850	! zoliqi = MAX(zoliq(i)zfice(i)-1.zice , 0.0)
	851	zoliqp(i) = MAX(zoliq(i)(1.-zfice(i))-1.zpluie , 0.0)
	852	zoliqi(i) = MAX(zoliq(i)zfice(i)-1.zice , 0.0)
[1472]	853	zoliq(i) = MAX(zoliq(i)-ztot , 0.0)
[1849]	854	!>AJ
[1472]	855	radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
	856	ENDIF
	857	ENDDO
	858	ENDDO
	859	!
[1849]	860	IF (.NOT. ice_thermo) THEN
	861	DO i = 1, klon
	862	IF (rneb(i,k).GT.0.0) THEN
[1472]	863	d_ql(i,k) = zoliq(i)
	864	zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) &
	865	* (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
[1849]	866	ENDIF
	867	ENDDO
	868	ELSE
	869	DO i = 1, klon
	870	IF (rneb(i,k).GT.0.0) THEN
	871	d_ql(i,k) = zoliq(i)
	872	!AJ<
	873	zrfl(i) = zrfl(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) &
	874	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
	875	zifl(i) = zifl(i)+ MAX(zcond(i)*zfice(i)-zoliqi(i),0.0) &
	876	(paprs(i,k)-paprs(i,k+1))/(RGdtime)
	877	! zrfl(i) = zrfl(i)+ zpluie &
	878	! (paprs(i,k)-paprs(i,k+1))/(RGdtime)
	879	! zifl(i) = zifl(i)+ zice &
	880	! (paprs(i,k)-paprs(i,k+1))/(RGdtime)
	881
	882	ENDIF
	883	ENDDO
	884	ENDIF
	885
	886	IF (ice_thermo) THEN
	887	DO i = 1, klon
	888	zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2
	889	zmelt = MIN(MAX(zmelt,0.),1.)
	890	zrfl(i)=zrfl(i)+zmelt*zifl(i)
	891	zifl(i)=zifl(i)*(1.-zmelt)
	892	! print*,zt(i),'octavio1'
	893	zt(i)=zt(i)-zifl(i)zmelt(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) &
	894	RLMLT/RCPD/(1.0+RVTMP2zq(i))
	895	! print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2'
	896	ENDDO
	897	ENDIF
	898
	899
	900	IF (.NOT. ice_thermo) THEN
	901	DO i = 1, klon
	902	IF (zt(i).LT.RTT) THEN
[1472]	903	psfl(i,k)=zrfl(i)
[1849]	904	ELSE
[1472]	905	prfl(i,k)=zrfl(i)
[1849]	906	ENDIF
	907	ENDDO
	908	ELSE
	909	! JAM*************************************************
	910	! Revoir partie ci-dessous: à quoi servent psfl et prfl?
	911	! *****************************************************
	912
	913	DO i = 1, klon
	914	! IF (zt(i).LT.RTT) THEN
	915	psfl(i,k)=zifl(i)
	916	! ELSE
	917	prfl(i,k)=zrfl(i)
	918	! ENDIF
	919	!>AJ
	920	ENDDO
	921	ENDIF
[1472]	922	!
[1849]	923	!
[1472]	924	! Calculer les tendances de q et de t:
	925	!
	926	DO i = 1, klon
	927	d_q(i,k) = zq(i) - q(i,k)
	928	d_t(i,k) = zt(i) - t(i,k)
	929	ENDDO
	930	!
	931	!AA--------------- Calcul du lessivage stratiforme -------------
[524]	932
[1472]	933	DO i = 1,klon
	934	!
[1742]	935	if(zcond(i).gt.zoliq(i)+1.e-10) then
	936	beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime
	937	else
	938	beta(i,k) = 0.
	939	endif
[1472]	940	zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) &
	941	* (paprs(i,k)-paprs(i,k+1))/RG
	942	IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
	943	!AA lessivage nucleation LMD5 dans la couche elle-meme
[2006]	944	IF (iflag_t_glace.EQ.0) THEN
	945	if (t(i,k) .GE. t_glace_min_old) THEN
[1472]	946	zalpha_tr = a_tr_sca(3)
	947	else
	948	zalpha_tr = a_tr_sca(4)
	949	endif
[2006]	950	ELSE ! of IF (iflag_t_glace.EQ.0)
	951	if (t(i,k) .GE. t_glace_min) THEN
	952	zalpha_tr = a_tr_sca(3)
	953	else
	954	zalpha_tr = a_tr_sca(4)
	955	endif
	956	ENDIF
[1472]	957	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
	958	pfrac_nucl(i,k)=pfrac_nucl(i,k)(1.-zneb(i)zfrac_lessi)
	959	frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
	960	!
	961	! nucleation avec un facteur -1 au lieu de -0.5
	962	zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
	963	pfrac_1nucl(i,k)=pfrac_1nucl(i,k)(1.-zneb(i)zfrac_lessi)
	964	ENDIF
	965	!
	966	ENDDO ! boucle sur i
	967	!
	968	!AA Lessivage par impaction dans les couches en-dessous
	969	DO kk = k-1, 1, -1
[524]	970	DO i = 1, klon
[1472]	971	IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
[2006]	972	IF (iflag_t_glace.EQ.0) THEN
	973	if (t(i,kk) .GE. t_glace_min_old) THEN
[1472]	974	zalpha_tr = a_tr_sca(1)
	975	else
	976	zalpha_tr = a_tr_sca(2)
	977	endif
[2006]	978	ELSE ! of IF (iflag_t_glace.EQ.0)
	979	if (t(i,kk) .GE. t_glace_min) THEN
	980	zalpha_tr = a_tr_sca(1)
	981	else
	982	zalpha_tr = a_tr_sca(2)
	983	endif
	984	ENDIF
[1472]	985	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
	986	pfrac_impa(i,kk)=pfrac_impa(i,kk)(1.-zneb(i)zfrac_lessi)
	987	frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
	988	ENDIF
[524]	989	ENDDO
[1472]	990	ENDDO
	991	!
	992	!AA----------------------------------------------------------
	993	! FIN DE BOUCLE SUR K
	994	end DO
	995	!
	996	!AA-----------------------------------------------------------
	997	!
	998	! Pluie ou neige au sol selon la temperature de la 1ere couche
	999	!
	1000	DO i = 1, klon
	1001	IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
[1849]	1002	!AJ<
	1003	!! snow(i) = zrfl(i)
	1004	snow(i) = zrfl(i)+zifl(i)
	1005	!>AJ
[1472]	1006	zlh_solid(i) = RLSTT-RLVTT
	1007	ELSE
	1008	rain(i) = zrfl(i)
	1009	zlh_solid(i) = 0.
	1010	ENDIF
	1011	ENDDO
	1012	!
	1013	! For energy conservation : when snow is present, the solification
	1014	! latent heat is considered.
	1015	DO k = 1, klev
	1016	DO i = 1, klon
	1017	zcpair=RCPD(1.0+RVTMP2(q(i,k)+d_q(i,k)))
	1018	zmair=(paprs(i,k)-paprs(i,k+1))/RG
	1019	zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime
	1020	d_t(i,k) = d_t(i,k) + zlh_solid(i) zm_solid / (zcpairzmair)
	1021	END DO
	1022	END DO
	1023	!
[883]	1024
[1472]	1025	if (ncoreczq>0) then
[1575]	1026	WRITE(lunout,*)'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.'
[1472]	1027	endif
	1028
	1029	END SUBROUTINE fisrtilp

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

Download in other formats:

Original Format