Context Navigation

fisrtilp.F90 @ 1893

Last change on this file since 1893 was 1864, checked in by Laurent Fairhead, 11 years ago

Création d'une nouvelle testing:

merge des modifications du trunk entre r1796 et r1860

New testing version

merged modifications between r1796 and r1860 from the trunk

i.e.
svn merge -r1796:1860 http://svn.lmd.jussieu.fr/LMDZ/LMDZ5/trunk

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 1864 2013-09-11 09:45:01Z fairhead $
[524]	3	!
[1472]	4	!
	5	SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs, &
[1750]	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, &
[1864]	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)
[1864]	54	REAL ztfondue, qsl, qsi
[1403]	55
[1472]	56	logical lognormale(klon)
[1864]	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
[1864]	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
[1864]	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)
[1864]	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
[1750]	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
[1864]	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	!
[1664]	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
[1664]	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.
[1750]	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
[1864]	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
[1864]	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
[1864]	292	! Calculer l'evaporation de la precipitation
	293	!
	294
	295
[1472]	296	IF (evap_prec) THEN
	297	DO i = 1, klon
[1864]	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
[1864]	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
[1864]	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
[1864]	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
[1864]	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
[1864]	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
[1796]	568	zcond(i) = MAX(0.0,zqn(i)-zqs(i))/(1+zdqs(i))
[1146]	569	rhcl(i,k)=1.0
	570	ELSE
[1796]	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
[1864]	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)
[1864]	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
[1864]	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)
[1864]	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)
[1864]	657	! Initialization of zpluie and zice:
	658	zpluie=0
	659	zice=0
[1472]	660	IF (zneb(i).EQ.seuil_neb) THEN
	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))
[1864]	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))
[1864]	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)
[1864]	696	!>AJ
[1472]	697	radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1)
	698	ENDIF
	699	ENDDO
	700	ENDDO
	701	!
[1864]	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)
[1864]	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)
[1864]	746	ELSE
[1472]	747	prfl(i,k)=zrfl(i)
[1864]	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	!
[1864]	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	!
[1750]	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
[1864]	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
[1664]	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