Context Navigation

fisrtilp.F90 @ 1849

Last change on this file since 1849 was 1849, checked in by Ehouarn Millour, 11 years ago
Including the thermodynamics of ice in the convection scheme (iactive only if iflag_ice_thermo==1). CR+JYG
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 26.5 KB

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

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

Download in other formats:

Original Format