Context Navigation

fisrtilp_tr.F @ 293

Last change on this file since 293 was 230, checked in by lmdzadmin, 23 years ago
Merge de la physique avec la branche principale LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 12.3 KB

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

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

Download in other formats:

Original Format