Context Navigation

fisrtilp_tr.F @ 191

Last change on this file since 191 was 60, checked in by lmdz, 25 years ago
Calcul de prfl et psfl en ciel clair aussi. O. Boucher LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 12.5 KB

Rev	Line
[2]	1	SUBROUTINE fisrtilp_tr(dtime,paprs,pplay,t,q,
	2	s d_t, d_q, d_ql, rneb, radliq, rain, snow,
	3	s pfrac_impa, pfrac_nucl, pfrac_1nucl,
[27]	4	s frac_impa, frac_nucl,
	5	s prfl, psfl)
[2]	6
	7	c
	8	IMPLICIT none
	9	c======================================================================
	10	c Auteur(s): Z.X. Li (LMD/CNRS)
	11	c Date: le 20 mars 1995
	12	c Objet: condensation et precipitation stratiforme.
	13	c schema de nuage
	14	c======================================================================
	15	c======================================================================
	16	#include "dimensions.h"
	17	#include "dimphy.h"
	18	#include "YOMCST.h"
	19	#include "tracstoke.h"
	20	c
	21	c Arguments:
	22	c
	23	REAL dtime ! intervalle du temps (s)
	24	REAL paprs(klon,klev+1) ! pression a inter-couche
	25	REAL pplay(klon,klev) ! pression au milieu de couche
	26	REAL t(klon,klev) ! temperature (K)
	27	REAL q(klon,klev) ! humidite specifique (kg/kg)
	28	REAL d_t(klon,klev) ! incrementation de la temperature (K)
	29	REAL d_q(klon,klev) ! incrementation de la vapeur d'eau
	30	REAL d_ql(klon,klev) ! incrementation de l'eau liquide
	31	REAL rneb(klon,klev) ! fraction nuageuse
	32	REAL radliq(klon,klev) ! eau liquide utilisee dans rayonnements
	33	REAL rain(klon) ! pluies (mm/s)
	34	REAL snow(klon) ! neige (mm/s)
[27]	35	REAL prfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
	36	REAL psfl(klon,klev+1) ! flux d'eau precipitante aux interfaces (kg/m2/s)
[2]	37	cAA
	38	c Coeffients de fraction lessivee : pour OFF-LINE
	39	c
	40	REAL pfrac_nucl(klon,klev)
	41	REAL pfrac_1nucl(klon,klev)
	42	REAL pfrac_impa(klon,klev)
	43	c
	44	c Fraction d'aerosols lessivee par impaction et par nucleation
	45	c POur ON-LINE
	46	c
	47	REAL frac_impa(klon,klev)
	48	REAL frac_nucl(klon,klev)
	49	cAA
	50	c
	51	c Options du programme:
	52	c
	53	REAL seuil_neb ! un nuage existe vraiment au-dela
	54	PARAMETER (seuil_neb=0.001)
	55	REAL ct ! inverse du temps pour qu'un nuage precipite
	56	PARAMETER (ct=1./1800.)
	57	REAL cl ! seuil de precipitation
[43]	58	PARAMETER (cl=2.6e-4)
	59	ccc PARAMETER (cl=2.3e-4)
	60	ccc PARAMETER (cl=2.0e-4)
[2]	61	INTEGER ninter ! sous-intervals pour la precipitation
	62	PARAMETER (ninter=5)
	63	LOGICAL evap_prec ! evaporation de la pluie
	64	PARAMETER (evap_prec=.TRUE.)
	65	REAL coef_eva
	66	PARAMETER (coef_eva=2.0E-05)
	67	LOGICAL calcrat ! calculer ratqs au lieu de fixer sa valeur
	68	REAL ratqs ! determine la largeur de distribution de vapeur
[43]	69	PARAMETER (calcrat=.TRUE.)
[2]	70	REAL zx_min, rat_max
	71	PARAMETER (zx_min=1.0, rat_max=0.01)
	72	REAL zx_max, rat_min
	73	PARAMETER (zx_max=0.1, rat_min=0.3)
	74	REAL zx
	75	c
	76	LOGICAL cpartiel ! condensation partielle
	77	PARAMETER (cpartiel=.TRUE.)
	78	REAL t_coup
	79	PARAMETER (t_coup=234.0)
	80	c
	81	c Variables locales:
	82	c
	83	INTEGER i, k, n, kk
	84	REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5
	85	REAL zrfl(klon), zrfln(klon), zqev, zqevt
	86	REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq
	87	REAL ztglace, zt(klon)
	88	INTEGER nexpo ! exponentiel pour glace/eau
	89	REAL zdz(klon),zrho(klon),ztot(klon), zrhol(klon)
	90	REAL zchau(klon),zfroi(klon),zfice(klon),zneb(klon)
	91	c
	92	LOGICAL appel1er
	93	SAVE appel1er
	94	c
	95	c---------------------------------------------------------------
	96	c
	97	cAA Variables traceurs:
	98	cAA Provisoire !!! Parametres alpha du lessivage
	99	cAA A priori on a 4 scavenging # possibles
	100	c
	101	REAL a_tr_sca(4)
	102	save a_tr_sca
	103	c
	104	c Variables intermediaires
	105	c
	106	REAL zalpha_tr
	107	REAL zfrac_lessi
[16]	108	REAL zprec_cond(klon)
[2]	109	cAA
	110	c---------------------------------------------------------------
	111	c
	112	c Fonctions en ligne:
	113	c
	114	REAL fallv ! vitesse de chute pour crystaux de glace
	115	REAL zzz
	116	#include "YOETHF.h"
	117	#include "FCTTRE.h"
[43]	118	fallv (zzz) = 3.29/2.0 * ((zzz)**0.16)
	119	ccc fallv (zzz) = 3.29/3.0 * ((zzz)**0.16)
	120	ccc fallv (zzz) = 3.29 * ((zzz)**0.16)
[2]	121	c
	122	DATA appel1er /.TRUE./
	123	c
	124	IF (appel1er) THEN
	125	c
	126	PRINT*, 'fisrtilp, calcrat:', calcrat
	127	PRINT*, 'fisrtilp, ninter:', ninter
	128	PRINT*, 'fisrtilp, evap_prec:', evap_prec
	129	PRINT*, 'fisrtilp, cpartiel:', cpartiel
	130	IF (ABS(dtime/FLOAT(ninter)-360.0).GT.0.001) THEN
	131	PRINT*, 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime
	132	PRINT*, 'Je prefere un sous-intervalle de 6 minutes'
	133	CALL abort
	134	ENDIF
	135	appel1er = .FALSE.
	136	c
	137	cAA initialiation provisoire
	138	a_tr_sca(1) = -0.5
	139	a_tr_sca(2) = -0.5
	140	a_tr_sca(3) = -0.5
	141	a_tr_sca(4) = -0.5
	142	c
	143	cAA Initialisation a 1 des coefs des fractions lessivees
	144	c
	145	DO k = 1, klev
	146	DO i = 1, klon
	147	pfrac_nucl(i,k)=1.
	148	pfrac_1nucl(i,k)=1.
	149	pfrac_impa(i,k)=1.
	150	ENDDO
	151	ENDDO
	152
	153	ENDIF ! test sur appel1er
	154	c
	155	cMAf Initialisation a 0 de zoliq
	156	DO i = 1, klon
	157	zoliq(i)=0.
	158	ENDDO
	159	c Determiner les nuages froids par leur temperature
	160	c
	161	ztglace = RTT - 15.0
	162	nexpo = 6
	163	ccc nexpo = 1
	164	c
	165	c Initialiser les sorties:
	166	c
[27]	167	DO k = 1, klev+1
	168	DO i = 1, klon
	169	prfl(i,k) = 0.0
	170	psfl(i,k) = 0.0
	171	ENDDO
	172	ENDDO
	173
[2]	174	DO k = 1, klev
	175	DO i = 1, klon
	176	d_t(i,k) = 0.0
	177	d_q(i,k) = 0.0
	178	d_ql(i,k) = 0.0
	179	rneb(i,k) = 0.0
	180	radliq(i,k) = 0.0
	181	frac_nucl(i,k) = 1.
	182	frac_impa(i,k) = 1.
	183	ENDDO
	184	ENDDO
	185	DO i = 1, klon
	186	rain(i) = 0.0
	187	snow(i) = 0.0
	188	ENDDO
	189	c
	190	c Initialiser le flux de precipitation a zero
	191	c
	192	DO i = 1, klon
	193	zrfl(i) = 0.0
	194	zneb(i) = seuil_neb
	195	ENDDO
	196	c
	197	c
	198	cAA Pour plus de securite
	199
	200	zalpha_tr = 0.
	201	zfrac_lessi = 0.
	202
	203	cAA----------------------------------------------------------
	204	c
	205	c Boucle verticale (du haut vers le bas)
	206	c
	207	DO 9999 k = klev, 1, -1
	208	c
	209	cAA----------------------------------------------------------
	210	c
	211	DO i = 1, klon
	212	zt(i)=t(i,k)
	213	zq(i)=q(i,k)
	214	ENDDO
	215	c
	216	c Calculer l'evaporation de la precipitation
	217	c
	218	IF (evap_prec) THEN
	219	DO i = 1, klon
	220	IF (zrfl(i) .GT.0.) THEN
	221	IF (thermcep) THEN
	222	zdelta=MAX(0.,SIGN(1.,RTT-zt(i)))
	223	zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
	224	zqs(i)=MIN(0.5,zqs(i))
	225	zcor=1./(1.-RETV*zqs(i))
	226	zqs(i)=zqs(i)*zcor
	227	ELSE
	228	IF (zt(i) .LT. t_coup) THEN
	229	zqs(i) = qsats(zt(i)) / pplay(i,k)
	230	ELSE
	231	zqs(i) = qsatl(zt(i)) / pplay(i,k)
	232	ENDIF
	233	ENDIF
	234	zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) )
	235	zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i))
	236	. * (paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)RD/RG
	237	zqevt = MAX(0.0,MIN(zqevt,zrfl(i)))
	238	. * RG*dtime/(paprs(i,k)-paprs(i,k+1))
	239	zqev = MIN (zqev, zqevt)
	240	zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1))
	241	. /RG/dtime
	242	zq(i) = zq(i) - (zrfln(i)-zrfl(i))
	243	. * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
	244	zt(i) = zt(i) + (zrfln(i)-zrfl(i))
	245	. * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
	246	. * RLVTT/RCPD/(1.0+RVTMP2*zq(i))
	247	zrfl(i) = zrfln(i)
	248	ENDIF
	249	ENDDO
	250	ENDIF
	251	c
	252	c Calculer Qs et L/Cp*dQs/dT:
	253	c
	254	IF (thermcep) THEN
	255	DO i = 1, klon
	256	zdelta = MAX(0.,SIGN(1.,RTT-zt(i)))
	257	zcvm5 = R5LESRLVTT(1.-zdelta) + R5IESRLSTTzdelta
	258	zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i))
	259	zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k)
	260	zqs(i) = MIN(0.5,zqs(i))
	261	zcor = 1./(1.-RETV*zqs(i))
	262	zqs(i) = zqs(i)*zcor
	263	zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor)
	264	ENDDO
	265	ELSE
	266	DO i = 1, klon
	267	IF (zt(i).LT.t_coup) THEN
	268	zqs(i) = qsats(zt(i))/pplay(i,k)
	269	zdqs(i) = dqsats(zt(i),zqs(i))
	270	ELSE
	271	zqs(i) = qsatl(zt(i))/pplay(i,k)
	272	zdqs(i) = dqsatl(zt(i),zqs(i))
	273	ENDIF
	274	ENDDO
	275	ENDIF
	276	c
	277	c Determiner la condensation partielle et calculer la quantite
	278	c de l'eau condensee:
	279	c
	280	IF (cpartiel) THEN
	281	DO i = 1, klon
	282	c
	283	zx = pplay(i,k)/paprs(i,1)
	284	zx = (zx_max-zx)/(zx_max-zx_min)
	285	zx = MIN(MAX(zx,0.0),1.0)
	286	zx = zx * zx * zx
	287	ratqs = zx * (rat_max-rat_min) + rat_min
[43]	288	IF (.NOT.calcrat) ratqs=0.05
[2]	289	c
	290	zdelq = ratqs * zq(i)
	291	rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq)
	292	zqn(i) = (zq(i)+zdelq+zqs(i))/2.0
	293	IF (rneb(i,k) .LE. 0.0) zqn(i) = 0.0
	294	IF (rneb(i,k) .GE. 1.0) zqn(i) = zq(i)
	295	rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
	296	zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
	297	ENDDO
	298	ELSE
	299	DO i = 1, klon
	300	IF (zq(i).GT.zqs(i)) THEN
	301	rneb(i,k) = 1.0
	302	ELSE
	303	rneb(i,k) = 0.0
	304	ENDIF
	305	zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i))
	306	ENDDO
	307	ENDIF
	308	c
	309	DO i = 1, klon
	310	zq(i) = zq(i) - zcond(i)
	311	zt(i) = zt(i) + zcond(i) * RLVTT/RCPD
	312	ENDDO
	313	c
	314	c Partager l'eau condensee en precipitation et eau liquide nuageuse
	315	c
	316	DO i = 1, klon
	317	IF (rneb(i,k).GT.0.0) THEN
	318	zoliq(i) = zcond(i)
	319	zrho(i) = pplay(i,k) / zt(i) / RD
	320	zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG)
	321	zfice(i) = 1.0 - (zt(i)-ztglace) / (273.13-ztglace)
	322	zfice(i) = MIN(MAX(zfice(i),0.0),1.0)
	323	zfice(i) = zfice(i)**nexpo
	324	zneb(i) = MAX(rneb(i,k), seuil_neb)
	325	radliq(i,k) = zoliq(i)/FLOAT(ninter+1)
	326	ENDIF
	327	ENDDO
	328	c
	329	DO n = 1, ninter
	330	DO i = 1, klon
	331	IF (rneb(i,k).GT.0.0) THEN
	332	zchau(i) = ctdtime/FLOAT(ninter) zoliq(i)
	333	. * (1.0-EXP(-(zoliq(i)/zneb(i)/cl)*2)) (1.-zfice(i))
	334	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
	335	zfroi(i) = dtime/FLOAT(ninter)/zdz(i)*zoliq(i)
	336	. fallv(zrhol(i)) zfice(i)
	337	ztot(i) = zchau(i) + zfroi(i)
	338	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
	339	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
	340	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
	341	radliq(i,k) = radliq(i,k) + zoliq(i)/FLOAT(ninter+1)
	342	ENDIF
	343	ENDDO
	344	ENDDO
	345	c
	346	DO i = 1, klon
	347	IF (rneb(i,k).GT.0.0) THEN
	348	d_ql(i,k) = zoliq(i)
	349	zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0)
	350	. * (paprs(i,k)-paprs(i,k+1))/(RG*dtime)
	351	ENDIF
[60]	352	IF (zt(i).LT.RTT) THEN
	353	psfl(i,k)=zrfl(i)
	354	ELSE
	355	prfl(i,k)=zrfl(i)
	356	ENDIF
[2]	357	ENDDO
	358	c
	359	c Calculer les tendances de q et de t:
	360	c
	361	DO i = 1, klon
	362	d_q(i,k) = zq(i) - q(i,k)
	363	d_t(i,k) = zt(i) - t(i,k)
	364	ENDDO
	365	c
	366	cAA--------------- Calcul du lessivage stratiforme -------------
	367
	368	DO i = 1,klon
	369	c
[16]	370	zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)
[2]	371	. * (paprs(i,k)-paprs(i,k+1))/RG
[16]	372	IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
[2]	373	cAA lessivage nucleation LMD5 dans la couche elle-meme
	374	if (t(i,k) .GE. ztglace) THEN
	375	zalpha_tr = a_tr_sca(3)
	376	else
	377	zalpha_tr = a_tr_sca(4)
	378	endif
[16]	379	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
[2]	380	pfrac_nucl(i,k)=pfrac_nucl(i,k)(1.-zneb(i)zfrac_lessi)
	381	frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi
	382	c
	383	c nucleation avec un facteur -1 au lieu de -0.5
[16]	384	zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i))
[2]	385	pfrac_1nucl(i,k)=pfrac_1nucl(i,k)(1.-zneb(i)zfrac_lessi)
	386	ENDIF
	387	c
	388	ENDDO ! boucle sur i
	389	c
[16]	390	cAA Lessivage par impaction dans les couches en-dessous
	391	DO kk = k-1, 1, -1
	392	DO i = 1, klon
	393	IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN
	394	if (t(i,kk) .GE. ztglace) THEN
	395	zalpha_tr = a_tr_sca(1)
	396	else
	397	zalpha_tr = a_tr_sca(2)
	398	endif
	399	zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i))
	400	pfrac_impa(i,kk)=pfrac_impa(i,kk)(1.-zneb(i)zfrac_lessi)
	401	frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi
	402	ENDIF
	403	ENDDO
	404	ENDDO
	405	c
[2]	406	cAA----------------------------------------------------------
	407	c FIN DE BOUCLE SUR K
	408	9999 CONTINUE
	409	c
	410	cAA-----------------------------------------------------------
	411	c
	412	c Pluie ou neige au sol selon la temperature de la 1ere couche
	413	c
	414	DO i = 1, klon
	415	IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN
	416	snow(i) = zrfl(i)
	417	ELSE
	418	rain(i) = zrfl(i)
	419	ENDIF
	420	ENDDO
	421	c
	422	RETURN
	423	END

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

Download in other formats:

Original Format