Context Navigation

cvltr_scav.F90 @ 5403

Last change on this file since 5403 was 2408, checked in by Laurent Fairhead, 9 years ago
Merged trunk changes r2298:2396 into testing branch
File size: 28.8 KB

Rev	Line
[2147]	1	!
	2	! $Id $
	3	!
	4	SUBROUTINE cvltr_scav(pdtime, da, phi,phi2,d1a,dam, mpIN,epIN, &
	5	sigd,sij,wght_cvfd,clw,elij,epmlmMm,eplaMm, &
	6	pmflxrIN,pmflxsIN,ev,te,wdtrainA,wdtrainM, &
	7	paprs,it,tr,upd,dnd,inb,icb, &
	8	ccntrAA_3d,ccntrENV_3d,coefcoli_3d, &
	9	dtrcv,trsptd,dtrSscav,dtrsat,dtrUscav,qDi,qPr, &
	10	qPa,qMel,qTrdi,dtrcvMA,Mint, &
	11	zmfd1a,zmfphi2,zmfdam)
	12	!
	13	USE IOIPSL
	14	USE dimphy
[2408]	15	USE infotrac_phy, ONLY : nbtr,tname
[2147]	16	IMPLICIT NONE
	17	!=====================================================================
	18	! Objet : convection des traceurs / KE
	19	! Auteurs: M-A Filiberti and J-Y Grandpeix
	20	! modifiee par R Pilon : lessivage des traceurs / KE
	21	!=====================================================================
	22
	23	include "YOMCST.h"
	24	include "YOECUMF.h"
	25	include "conema3.h"
	26	include "chem.h"
	27
	28	! Entree
	29	REAL,INTENT(IN) :: pdtime
	30	REAL,DIMENSION(klon,klev),INTENT(IN) :: da
	31	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi
	32	! RomP
	33	REAL,DIMENSION(klon,klev),INTENT(IN) :: d1a,dam ! matrices pour simplifier
	34	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi2 ! l'ecriture des tendances
	35	!
	36	REAL,DIMENSION(klon,klev),INTENT(IN) :: mpIN
	37	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression aux 1/2 couches (bas en haut)
	38	INTEGER,INTENT(IN) :: it ! numero du traceur
	39	REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr ! q de traceur (bas en haut)
	40	REAL,DIMENSION(klon,klev),INTENT(IN) :: upd ! saturated updraft mass flux
	41	REAL,DIMENSION(klon,klev),INTENT(IN) :: dnd ! saturated downdraft mass flux
	42	!
	43	REAL,DIMENSION(klon,klev),INTENT(IN) :: wdtrainA ! masses precipitantes de l'asc adiab
	44	REAL,DIMENSION(klon,klev),INTENT(IN) :: wdtrainM ! masses precipitantes des melanges
	45	!JE REAL,DIMENSION(klon,klev),INTENT(IN) :: pmflxrIN ! vprecip: eau
	46	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: pmflxrIN ! vprecip: eau
	47	!JE REAL,DIMENSION(klon,klev),INTENT(IN) :: pmflxsIN ! vprecip: neige
	48	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: pmflxsIN ! vprecip: neige
	49	REAL,DIMENSION(klon,klev),INTENT(IN) :: ev ! evaporation cv30_routine
	50	REAL,DIMENSION(klon,klev),INTENT(IN) :: epIN
	51	REAL,DIMENSION(klon,klev),INTENT(IN) :: te
	52	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: sij ! fraction dair de lenv
	53	REAL,DIMENSION(klon,klev),INTENT(IN) :: wght_cvfd ! weights of the layers feeding convection
	54	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: elij ! contenu en eau condensée spécifique/conc deau condensée massique
	55	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: epmlmMm ! eau condensee precipitee dans mel masse dair sat
	56	REAL,DIMENSION(klon,klev),INTENT(IN) :: eplaMm ! eau condensee precipitee dans aa masse dair sat
	57
	58	REAL,DIMENSION(klon,klev),INTENT(IN) :: clw ! contenu en eau condensée dans lasc adiab
	59	REAL,DIMENSION(klon),INTENT(IN) :: sigd
	60	INTEGER,DIMENSION(klon),INTENT(IN) :: icb,inb
	61	!
	62	REAL,DIMENSION(klon,klev),INTENT(IN) :: ccntrAA_3d
	63	REAL,DIMENSION(klon,klev),INTENT(IN) :: ccntrENV_3d
	64	REAL,DIMENSION(klon,klev),INTENT(IN) :: coefcoli_3d
	65	!
	66	! Sortie
	67	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrcv ! tendance totale (bas en haut)
	68	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrcvMA ! M-A Filiberti
	69	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: trsptd ! tendance du transport
	70	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrSscav ! tendance du lessivage courant sat
	71	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrsat ! tendance trsp+sat scav
	72	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrUscav ! tendance du lessivage courant unsat
	73	!
	74	! Variables locales
	75	INTEGER :: i,j,k
	76	REAL,DIMENSION(klon,klev) :: dxpres ! difference de pression entre niveau (j+1) et (j)
	77	REAL :: pdtimeRG ! pas de temps * gravite
	78	REAL,DIMENSION(klon,nbtr) :: qfeed ! tracer concentration feeding convection
	79	! variables pour les courants satures
	80	REAL,DIMENSION(klon,klev,klev) :: zmd
	81	REAL,DIMENSION(klon,klev,klev) :: za
	82	REAL,DIMENSION(klon,klev,nbtr) :: zmfd,zmfa
	83	REAL,DIMENSION(klon,klev,nbtr) :: zmfp,zmfu
	84
	85	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfd1a
	86	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfdam
	87	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfphi2
	88
	89	! RomP ! les variables sont nettoyees des valeurs aberrantes
	90	REAL,DIMENSION(klon,klev) :: Pa, Pm ! pluie AA et mélanges, var temporaire
	91	REAL,DIMENSION(klon,klev) :: pmflxs,pmflxr ! pmflxrIN,pmflxsIN sans valeur aberante
	92	REAL,DIMENSION(klon,klev) :: mp ! flux de masse
	93	REAL,DIMENSION(klon,klev) :: ep ! fraction d'eau convertie en precipitation
	94	REAL,DIMENSION(klon,klev) :: evap ! evaporation : variable temporaire
	95	REAL,DIMENSION(klon,klev) :: rho !environmental density
	96
	97	REAL,DIMENSION(klon,klev) :: kappa ! denominateur du au calcul de la matrice
	98	! pour obtenir qd et qp
	99	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qTrdi ! traceurs descente air insature transport MA
	100	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qDi ! traceurs descente insaturees
	101	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qPr ! traceurs colonne precipitante
	102	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qPa ! traceurs dans les precip issues lasc. adiab.
	103	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qMel ! traceurs dans les precip issues des melanges
	104	REAL,DIMENSION(klon,klev,nbtr) :: qMeltmp ! variable temporaire
	105	REAL,DIMENSION(klon,klev,nbtr) :: qpmMint
	106	REAL,DIMENSION(klon,klev),INTENT(OUT) :: Mint
	107	! tendances
	108	REAL :: tdcvMA ! terme de transport de traceur (schema Marie Angele)
	109	REAL :: trsptrac ! terme de transport de traceur par l'air
	110	REAL :: scavtrac ! terme de lessivage courant sature
	111	REAL :: uscavtrac ! terme de lessivage courant insature
	112	! impaction
	113	!!! Correction apres discussion Romain P. / Olivier B.
	114	!!! REAL,PARAMETER :: rdrop=2.5e-3 ! rayon des gouttes d'eau
	115	REAL,PARAMETER :: rdrop=1.e-3 ! rayon des gouttes d'eau
	116	!!!
	117	REAL,DIMENSION(klon,klev) :: imp ! coefficient d'impaction
	118	!
	119	LOGICAL,DIMENSION(klon,klev) :: NO_precip
	120	! var tmp tests
	121	REAL :: conserv
	122	real :: conservMA
	123
[2298]	124	!jyg<
	125	!! ! ======================================================
	126	!! ! calcul de l'impaction
	127	!! ! ======================================================
	128	!!
	129	!! ! impaction sur la surface de la colonne de la descente insaturee
	130	!! ! On prend la moyenne des precip entre le niveau i+1 et i
	131	!! ! I=3/4* (P(1+1)+P(i))/2 / (sigdrrho_l)
	132	!! ! 1000kg/m3= densite de l'eau
	133	!! ! 0.75e-3 = 3/4 /1000
	134	!! ! Par la suite, I est tout le temps multiplie par sig_d pour avoir l'impaction sur la surface de la maille
	135	!!!! ! on le neglige ici pour simplifier le code
	136	!!
	137	!! DO j=1,klev-1
	138	!! DO i=1,klon
	139	!! imp(i,j) = coefcoli_3d(i,j)0.75e-3/rdrop &
	140	!! 0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j))
	141	!! ENDDO
	142	!! ENDDO
	143	!>jyg
[2147]	144	!
	145	! initialisation pour flux de traceurs, td et autre
	146	!
	147	trsptrac = 0.
	148	scavtrac = 0.
	149	uscavtrac = 0.
	150	qfeed(:,it) = 0. !RL
	151	DO j=1,klev
	152	DO i=1,klon
	153	zmfd(i,j,it)=0.
	154	zmfa(i,j,it)=0.
	155	zmfu(i,j,it)=0.
	156	zmfp(i,j,it)=0.
	157	zmfphi2(i,j,it)=0.
	158	zmfd1a(i,j,it)=0.
	159	zmfdam(i,j,it)=0.
	160	qDi(i,j,it)=0.
	161	qPr(i,j,it)=0.
	162	qPa(i,j,it)=0.
	163	qMel(i,j,it)=0.
	164	qMeltmp(i,j,it)=0.
	165	qTrdi(i,j,it)=0.
	166	kappa(i,j)=0.
	167	trsptd(i,j,it)=0.
	168	dtrsat(i,j,it)=0.
	169	dtrSscav(i,j,it)=0.
	170	dtrUscav(i,j,it)=0.
	171	dtrcv(i,j,it)=0.
	172	dtrcvMA(i,j,it)=0.
	173	evap(i,j)=0.
	174	dxpres(i,j)=0.
	175	qpmMint(i,j,it)=0.
	176	Mint(i,j)=0.
	177	END DO
	178	END DO
	179
	180	! suppression des valeurs très faibles (~1e-320)
	181	! multiplication de levaporation pour lavoir par unite de temps
	182	! et par unite de surface de la maille
	183	! -> cv30_unsat : evap : masse evaporee/s/(m2 de la descente)
	184	DO j=1,klev
	185	DO i=1,klon
	186	IF(ev(i,j).lt.1.e-16) THEN
	187	evap(i,j)=0.
	188	ELSE
	189	evap(i,j)=ev(i,j)*sigd(i)
	190	ENDIF
	191	END DO
	192	END DO
	193
	194	DO j=1,klev
	195	DO i=1,klon
	196	IF(j.LT.klev) THEN
	197	IF(epIN(i,j).LT.1.e-32) THEN
	198	ep(i,j)=0.
	199	ELSE
	200	ep(i,j)=epIN(i,j)
	201	ENDIF
	202	ELSE
	203	ep(i,j)=epmax
	204	ENDIF
	205	IF(mpIN(i,j).LT.1.e-32) THEN
	206	mp(i,j)=0.
	207	ELSE
	208	mp(i,j)=mpIN(i,j)
	209	ENDIF
	210	IF(pmflxsIN(i,j).LT.1.e-32) THEN
	211	pmflxs(i,j)=0.
	212	ELSE
	213	pmflxs(i,j)=pmflxsIN(i,j)
	214	ENDIF
	215	IF(pmflxrIN(i,j).LT.1.e-32) THEN
	216	pmflxr(i,j)=0.
	217	ELSE
	218	pmflxr(i,j)=pmflxrIN(i,j)
	219	ENDIF
	220	IF(wdtrainA(i,j).LT.1.e-32) THEN
	221	Pa(i,j)=0.
	222	ELSE
	223	Pa(i,j)=wdtrainA(i,j)
	224	ENDIF
	225	IF(wdtrainM(i,j).LT.1.e-32) THEN
	226	Pm(i,j)=0.
	227	ELSE
	228	Pm(i,j)=wdtrainM(i,j)
	229	ENDIF
	230	END DO
	231	END DO
	232
	233	!==========================================
	234	DO j = klev-1,1,-1
	235	DO i = 1,klon
	236	NO_precip(i,j) = (pmflxr(i,j+1)+pmflxs(i,j+1)).LT.1.e-10&
	237	.AND.Pa(i,j).LT.1.e-10.AND.Pm(i,j).LT.1.e-10
	238	END DO
	239	END DO
	240
[2298]	241	!jyg<
	242	! ======================================================
	243	! calcul de l'impaction
	244	! ======================================================
	245
	246	! impaction sur la surface de la colonne de la descente insaturee
	247	! On prend la moyenne des precip entre le niveau i+1 et i
	248	! I=3/4* (P(1+1)+P(i))/2 / (sigdrrho_l)
	249	! 1000kg/m3= densite de l'eau
	250	! 0.75e-3 = 3/4 /1000
	251	! Par la suite, I est tout le temps multiplie par sig_d pour avoir l'impaction sur la surface de la maille
	252	! on le neglige ici pour simplifier le code
	253
	254	DO j=1,klev-1
	255	DO i=1,klon
	256	imp(i,j) = coefcoli_3d(i,j)0.75e-3/rdrop &
	257	0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j))
	258	ENDDO
	259	ENDDO
	260	!>jyg
[2147]	261	! =========================================
	262	! calcul des tendances liees au downdraft
	263	! =========================================
	264	!cdir collapse
	265	DO k=1,klev
	266	DO j=1,klev
	267	DO i=1,klon
	268	zmd(i,j,k)=0.
	269	za (i,j,k)=0.
	270	END DO
	271	END DO
	272	END DO
	273	! calcul de la matrice d echange
	274	! matrice de distribution de la masse entrainee en k
	275	! commmentaire RomP : mp > 0
	276	DO k=1,klev-1
	277	DO i=1,klon
	278	zmd(i,k,k)=max(0.,mp(i,k)-mp(i,k+1)) ! ~ mk(k)
	279	END DO
	280	END DO
	281	DO k=2,klev
	282	DO j=k-1,1,-1
	283	DO i=1,klon
	284	IF(mp(i,j+1).GT.1.e-10) THEN
	285	zmd(i,j,k)=zmd(i,j+1,k)min(1.,mp(i,j)/mp(i,j+1)) !det ~ mk(j)=mk(j+1)mp(i,j)/mp(i,j+1)
	286	ENDIF
	287	END DO
	288	END DO
	289	END DO
	290	DO k=1,klev
	291	DO j=1,klev-1
	292	DO i=1,klon
	293	za(i,j,k)=max(0.,zmd(i,j+1,k)-zmd(i,j,k))
	294	END DO
	295	END DO
	296	END DO
	297	!!!!! quantite de traceur dans la descente d'air insaturee : 4 juin 2012
	298	DO k=1,klev
	299	DO j=1,klev-1
	300	DO i=1,klon
	301	IF(mp(i,j+1).GT.1.e-10) THEN
	302	qTrdi(i,j+1,it)=qTrdi(i,j+1,it)+(zmd(i,j+1,k)/mp(i,j+1))*tr(i,k,it)
	303	ELSE
	304	qTrdi(i,j,it)=0.!tr(i,j,it)
	305	ENDIF
	306	ENDDO
	307	ENDDO
	308	ENDDO
	309	!!!!!
	310	!
	311	! rajout du terme lie a l ascendance induite
	312	!
	313	DO j=2,klev
	314	DO i=1,klon
	315	za(i,j,j-1)=za(i,j,j-1)+mp(i,j)
	316	END DO
	317	END DO
	318	!
	319	! tendance courants insatures ! sans lessivage ancien schema
	320	!
	321	DO k=1,klev
	322	DO j=1,klev
	323	DO i=1,klon
	324	zmfd(i,j,it)=zmfd(i,j,it)+za(i,j,k)*(tr(i,k,it)-tr(i,j,it))
	325	END DO
	326	END DO
	327	END DO
	328	!
	329	! =========================================
	330	! calcul des tendances liees aux courants satures j <-> z ; k <-> z'
	331	! =========================================
	332	!RL
	333	! Feeding concentrations
	334	DO j=1,klev
	335	DO i=1,klon
	336	qfeed(i,it)=qfeed(i,it)+wght_cvfd(i,j)*tr(i,j,it)
	337	END DO
	338	END DO
	339	!RL
	340	!
	341	DO j=1,klev
	342	DO i=1,klon
	343	!RL
	344	!! zmfa(i,j,it)=da(i,j)*(tr(i,1,it)-tr(i,j,it)) ! da
	345	zmfa(i,j,it)=da(i,j)*(qfeed(i,it)-tr(i,j,it)) ! da
	346	!RL
	347	END DO
	348	END DO
	349	!
	350	DO k=1,klev
	351	DO j=1,klev
	352	DO i=1,klon
	353	zmfp(i,j,it)=zmfp(i,j,it)+phi(i,j,k)*(tr(i,k,it)-tr(i,j,it)) ! phi
	354	END DO
	355	END DO
	356	END DO
	357	! RomP ajout des matrices liees au lessivage
	358	DO j=1,klev
	359	DO i=1,klon
	360	zmfd1a(i,j,it)=d1a(i,j)*tr(i,1,it) ! da1
	361	zmfdam(i,j,it)=dam(i,j)*tr(i,1,it) ! dam
	362	END DO
	363	END DO
	364	DO k=1,klev
	365	DO j=1,klev
	366	DO i=1,klon
	367	zmfphi2(i,j,it)=zmfphi2(i,j,it)+phi2(i,j,k)*tr(i,k,it) ! psi
	368	END DO
	369	END DO
	370	END DO
	371	DO j=1,klev-1
	372	DO i=1,klon
	373	zmfu(i,j,it)=max(0.,upd(i,j+1)+dnd(i,j+1))*(tr(i,j+1,it)-tr(i,j,it))
	374	END DO
	375	END DO
	376	DO j=2,klev
	377	DO i=1,klon
	378	zmfu(i,j,it)=zmfu(i,j,it)+min(0.,upd(i,j)+dnd(i,j))*(tr(i,j,it)-tr(i,j-1,it))
	379	END DO
	380	END DO
	381	! ===================================================
	382	! calcul des tendances liees aux courants insatures
	383	! ===================================================
	384	! pression
	385	DO k=1, klev
	386	DO i=1, klon
	387	dxpres(i,k)=paprs(i,k)-paprs(i,k+1)
	388	ENDDO
	389	ENDDO
	390	pdtimeRG=pdtime*RG
	391
	392	! q_pa et q_pm traceurs issues des courants satures se retrouvant dans les precipitations
	393	DO j=1,klev
	394	DO i=1,klon
	395	IF(j.GE.icb(i).AND.j.LE.inb(i)) THEN
	396	IF(clw(i,j).GT.1.e-16) THEN
	397	!JE qPa(i,j,it)=ccntrAA_coef*tr(i,1,it)/clw(i,j)
	398	qPa(i,j,it)=ccntrAA_3d(i,j)*tr(i,1,it)/clw(i,j)
	399	ELSE
	400	qPa(i,j,it)=0.
	401	ENDIF
	402	ENDIF
	403	END DO
	404	END DO
	405
	406	! calcul de q_pm en 2 parties :
	407	! 1) calcul de sa valeur pour un niveau z' donne
	408	! 2) integration sur la verticale sur z'
	409	DO j=1,klev
	410	DO k=1,j-1
	411	DO i=1,klon
	412	IF(k.GE.icb(i).AND.k.LE.inb(i).AND.&
	413	j.LE.inb(i)) THEN
	414	IF(elij(i,k,j).GT.1.e-16) THEN
	415	!JE qMeltmp(i,j,it)=((1-ep(i,k))ccntrAA_coeftr(i,1,it)&
	416	!JE *(1.-sij(i,k,j)) +ccntrENV_coef&
	417	!JE tr(i,k,it)sij(i,k,j)) / elij(i,k,j)
	418	qMeltmp(i,j,it)=((1-ep(i,k))ccntrAA_3d(i,k)tr(i,1,it)&
	419	*(1.-sij(i,k,j)) +ccntrENV_3d(i,k)&
	420	tr(i,k,it)sij(i,k,j)) / elij(i,k,j)
	421	ELSE
	422	qMeltmp(i,j,it)=0.
	423	ENDIF
	424	qpmMint(i,j,it)=qpmMint(i,j,it) + qMeltmp(i,j,it)*epmlmMm(i,j,k)
	425	Mint(i,j)=Mint(i,j) + epmlmMm(i,j,k)
	426	ENDIF ! end if dans nuage
	427	END DO
	428	END DO
	429	END DO
	430
	431	DO j=1,klev
	432	DO i=1,klon
	433	IF(Mint(i,j).GT.1.e-16) THEN
	434	qMel(i,j,it)=qpmMint(i,j,it)/Mint(i,j)
	435	ELSE
	436	qMel(i,j,it)=0.
	437	ENDIF
	438	END DO
	439	END DO
	440
	441	! calcul de q_d et q_p traceurs de la descente precipitante
	442	DO j=klev-1,1,-1
	443	DO i=1,klon
	444	IF(mp(i,j+1).GT.mp(i,j).AND.mp(i,j+1).GT.1.e-10) THEN ! detrainement
	445	kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
	446	(-mp(i,j+1)-imp(i,j)/RG*dxpres(i,j))&
	447	+ (imp(i,j)/RGdxpres(i,j))(evap(i,j)/RG*dxpres(i,j)))
	448
	449	ELSEIF(mp(i,j).GT.mp(i,j+1).AND.mp(i,j).GT.1.e-10) THEN! entrainement
	450	IF(j.eq.1) THEN
	451	kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
	452	(-mp(i,2)-imp(i,j)/RG*dxpres(i,j))&
	453	+ (imp(i,j)/RGdxpres(i,j))(evap(i,j)/RG*dxpres(i,j)))
	454	ELSE
	455	kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
	456	(-mp(i,j)-imp(i,j)/RG*dxpres(i,j))&
	457	+ (imp(i,j)/RGdxpres(i,j))(evap(i,j)/RG*dxpres(i,j)))
	458	ENDIF
	459	ELSE
	460	kappa(i,j)=1.
	461	ENDIF
	462	ENDDO
	463	ENDDO
	464
	465	DO j=klev-1,1,-1
	466	DO i=1,klon
	467	IF (abs(kappa(i,j)).LT.1.e-25) THEN !si denominateur nul (il peut y avoir des mp!=0)
	468	kappa(i,j)=1.
	469	IF(j.eq.1) THEN
	470	qDi(i,j,it)=qDi(i,j+1,it) !orig tr(i,j,it) ! mp(1)=0 donc tout vient de la couche supérieure
	471	ELSEIF(mp(i,j+1).GT.mp(i,j).AND.mp(i,j+1).GT.1.e-10) THEN
	472	qDi(i,j,it)=qDi(i,j+1,it)
	473	ELSEIF(mp(i,j).GT.mp(i,j+1).AND.mp(i,j).GT.1.e-10) THEN! entrainement
	474	qDi(i,j,it)=(-mp(i,j+1)(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)tr(i,j,it))/(-mp(i,j))
	475	ELSE ! si mp (i)=0 et mp(j+1)=0
	476	qDi(i,j,it)=tr(i,j,it) ! orig 0.
	477	ENDIF
	478
	479	IF(NO_precip(i,j)) THEN
	480	qPr(i,j,it)=0.
	481	ELSE
	482	qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
	483	Pa(i,j)qPa(i,j,it)+Pm(i,j)qMel(i,j,it)&
	484	+imp(i,j)/RGdxpres(i,j)qDi(i,j,it))/&
	485	(pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))
	486	ENDIF
	487	ELSE ! denominateur non nul
	488	kappa(i,j)=1./kappa(i,j)
	489	! calcul de qd et qp
	490	!!jyg (20130119) correction pour le sommet du nuage
	491	!! if(j.GE.inb(i)) THEN !au-dessus du nuage, sommet inclu
	492	if(j.GT.inb(i)) THEN !au-dessus du nuage
	493	qDi(i,j,it)=tr(i,j,it) ! pas de descente => environnement = descente insaturee
	494	qPr(i,j,it)=0.
	495
	496	! vvv premiere couche du modele ou mp(1)=0 ! det tout le temps vvv
	497	ELSEIF(j.eq.1) THEN
	498	if(mp(i,2).GT.1.e-10) THEN !mp(2) non nul -> detrainement (car mp(1) = 0) !ent pas possible
	499	if(NO_precip(i,j)) THEN !pas de precip en (i)
	500	qDi(i,j,it)=qDi(i,j+1,it)
	501	qPr(i,j,it)=0.
	502	ELSE
	503	qDi(i,j,it)=kappa(i,j)*(&
	504	(-evap(i,j)/RGdxpres(i,j))((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
	505	Pa(i,j)qPa(i,j,it)+Pm(i,j)qMel(i,j,it)) +&
	506	(pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
	507	(-mp(i,j+1)*qDi(i,j+1,it)))
	508
	509	qPr(i,j,it)=kappa(i,j)*(&
	510	(-mp(i,j+1)-imp(i,j)/RGdxpres(i,j))&
	511	((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
	512	Pa(i,j)qPa(i,j,it)+Pm(i,j)qMel(i,j,it))&
	513	+(-mp(i,j+1)qDi(i,j+1,it)) (imp(i,j)/RG*dxpres(i,j)))
	514	ENDIF
	515
	516	ELSE !mp(2) nul -> plus de descente insaturee -> pluie agit sur environnement
	517	qDi(i,j,it)=tr(i,j,it) ! orig 0.
	518	if(NO_precip(i,j)) THEN
	519	qPr(i,j,it)=0.
	520	ELSE
	521	qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
	522	Pa(i,j)qPa(i,j,it)+Pm(i,j)qMel(i,j,it)&
	523	+imp(i,j)/RGdxpres(i,j)tr(i,j,it))/&
	524	(pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))
	525	ENDIF
	526
	527	ENDIF !mp(2) nul ou non
	528
	529	! vvv (j!=1.AND.j.LT.inb(i)) en-dessous du sommet nuage vvv
	530	ELSE
	531	!------------------------------------------------------------- detrainement
	532	if(mp(i,j+1).GT.mp(i,j).AND.mp(i,j+1).GT.1.e-10) THEN !mp(i,j).GT.1.e-10) THEN
	533	if(NO_precip(i,j)) THEN
	534	qDi(i,j,it)=qDi(i,j+1,it)
	535	qPr(i,j,it)=0.
	536	ELSE
	537	qDi(i,j,it)=kappa(i,j)*(&
	538	(-evap(i,j)/RGdxpres(i,j))((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
	539	Pa(i,j)qPa(i,j,it)+Pm(i,j)qMel(i,j,it)) +&
	540	(pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
	541	(-mp(i,j+1)*qDi(i,j+1,it)))
	542	!
	543	qPr(i,j,it)=kappa(i,j)*(&
	544	(-mp(i,j+1)-imp(i,j)/RGdxpres(i,j))&
	545	((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
	546	Pa(i,j)qPa(i,j,it)+Pm(i,j)qMel(i,j,it))&
	547	+(-mp(i,j+1)qDi(i,j+1,it)) (imp(i,j)/RG*dxpres(i,j)))
	548	ENDIF !precip
	549	!------------------------------------------------------------- entrainement
	550	ELSEIF(mp(i,j).GT.mp(i,j+1).AND.mp(i,j).GT.1.e-10) THEN
	551	if(NO_precip(i,j)) THEN
	552	qDi(i,j,it)=(-mp(i,j+1)(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)tr(i,j,it))/(-mp(i,j))
	553	qPr(i,j,it)=0.
	554	ELSE
	555	qDi(i,j,it)=kappa(i,j)*(&
	556	(-evap(i,j)/RGdxpres(i,j))((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
	557	Pa(i,j)qPa(i,j,it)+Pm(i,j)qMel(i,j,it)) +&
	558	(pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
	559	(-mp(i,j+1)(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)tr(i,j,it)))
	560	!
	561	qPr(i,j,it)=kappa(i,j)*(&
	562	(-mp(i,j)-imp(i,j)/RGdxpres(i,j))&
	563	((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
	564	Pa(i,j)qPa(i,j,it)+Pm(i,j)qMel(i,j,it))&
	565	+(-mp(i,j+1)(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)tr(i,j,it))*&
	566	(imp(i,j)/RG*dxpres(i,j)))
	567	ENDIF !precip
	568	!------------------------------------------------------------- ENDIF ! ent/det
	569	ELSE !mp nul
	570	qDi(i,j,it)=tr(i,j,it) ! orig 0.
	571	if(NO_precip(i,j)) THEN
	572	qPr(i,j,it)=0.
	573	ELSE
	574	qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
	575	Pa(i,j)qPa(i,j,it)+Pm(i,j)qMel(i,j,it)&
	576	+imp(i,j)/RGdxpres(i,j)tr(i,j,it))/&
	577	(pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))
	578	ENDIF
	579	ENDIF ! mp nul ou non
	580	ENDIF ! condition sur j
	581	ENDIF ! kappa
	582	ENDDO
	583	ENDDO
	584
	585	!! print test descente insaturee
	586	! DO j=klev,1,-1
	587	! DO i=1,klon
	588	! if(it.eq.3) THEN
	589	! write(*,'(I2,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12)') j,&
	590	!! 'zmfdam',zmfdam(i,j,it),'zmfpsi',zmfphi2(i,j,it),&
	591	! 'zmfdam+zmfpsi',zmfdam(i,j,it)+zmfphi2(i,j,it),'qpmMint',qpmMint(i,j,it),&
	592	! 'Pm',Pm(i,j),'Mint',Mint(i,j),&
	593	!! 'zmfa',zmfa(i,j,it),'zmfp',zmfp(i,j,it),&
	594	! 'zmfdam',zmfdam(i,j,it),'zmfpsi',zmfphi2(i,j,it),'zmfd1a',zmfd1a(i,j,it)
	595	!! 'Pa',Pa(i,j),'eplaMm',eplaMm(i,j)
	596	!! 'zmfd1a=da1qa',zmfd1a(i,j,it),'PaqPa',wdtrainA(i,j)*qPa(i,j,it),'da1',d1a(i,j)
	597	! ENDIF
	598	! ENDDO
	599	! ENDDO
	600
	601
	602	! ===================================================
	603	! calcul final des tendances
	604	! ===================================================
	605
	606	DO k=klev-1,1,-1
	607	DO i=1, klon
	608	! transport
	609	tdcvMA=zmfd(i,k,it)+zmfu(i,k,it)+zmfa(i,k,it)+zmfp(i,k,it) ! double comptage des downdraft insatures
	610	trsptrac=zmfu(i,k,it)+zmfa(i,k,it)+zmfp(i,k,it)
	611	! lessivage courants satures
	612	!JE scavtrac=-ccntrAA_coef*zmfd1a(i,k,it)&
	613	!JE -zmfphi2(i,k,it)*ccntrENV_coef&
	614	!JE -zmfdam(i,k,it)*ccntrAA_coef
	615	scavtrac=-ccntrAA_3d(i,k)*zmfd1a(i,k,it)&
	616	-zmfphi2(i,k,it)*ccntrENV_3d(i,k)&
	617	-zmfdam(i,k,it)*ccntrAA_3d(i,k)
	618	! lessivage courants insatures
	619	if(k.LE.inb(i).AND.k.GT.1) THEN ! tendances dans le nuage
	620	!------------------------------------------------------------- detrainement
	621	if(mp(i,k+1).GT.mp(i,k).AND.mp(i,k+1).GT.1.e-10) THEN
	622	uscavtrac= (-mp(i,k)+mp(i,k+1))*(qDi(i,k,it)-tr(i,k,it))&
	623	+ mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))
	624	!
	625	! if(it.eq.3) write(*,'(I2,1X,a,5X,e20.12,82X,a,e20.12)')k,' det incloud',&
	626	! (-mp(i,k)+mp(i,k+1))(qDi(i,k,it)-tr(i,k,it))pdtimeRG/dxpres(i,k)+&
	627	! mp(i,k)(tr(i,k-1,it)-tr(i,k,it))pdtimeRG/dxpres(i,k),&
	628	! 'mp',mp(i,k)
	629	!------------------------------------------------------------- entrainement
	630	ELSEIF(mp(i,k).GT.mp(i,k+1).AND.mp(i,k).GT.1.e-10) THEN
	631	uscavtrac= mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))
	632	!
	633	! if(it.eq.3) write(,'(I2,1X,a,5X,e20.12,82X,a,e20.12)')k,' ent incloud',uscavtracpdtimeRG/dxpres(i,k), 'mp',mp(i,k)
	634	!=!------------------------------------------------------------- end ent/det
	635	ELSE ! mp(i,k+1)=0. et mp(i,k)=0. pluie directement sur l environnement
	636
	637	if(NO_precip(i,k)) THEN
	638	uscavtrac=0.
	639	! if(it.eq.3) write(,'(I2,1X,a,e20.12,82X,a,e20.12)')k,' no P ent incloud',uscavtracpdtimeRG/dxpres(i,k), 'mp',mp(i,k)
	640	ELSE
	641	uscavtrac=-imp(i,k)tr(i,k,it)dxpres(i,k)/RG+evap(i,k)qPr(i,k,it)dxpres(i,k)/RG
	642	! if(it.eq.3) write(,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtracpdtimeRG/dxpres(i,k), 'mp',mp(i,k)
	643	!!JE adds
	644	! if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'imp',imp(i,k)
	645	! if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'tr',tr(i,k,it)
	646	! if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'evap',evap(i,k)
	647	! if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'qPr',qPr(i,k,it)
	648	! if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'dxpres',dxpres(i,k)
	649	!!Je end
	650
	651	ENDIF
	652	ENDIF ! mp/det/ent
	653	!------------------------------------------------------------- premiere couche
	654	ELSEIF(k.eq.1) THEN ! mp(1)=0.
	655	if(mp(i,2).GT.1.e-10) THEN !detrainement
	656	uscavtrac= (-0.+mp(i,2))*(qDi(i,k,it)-tr(i,k,it)) !&
	657	! + mp(i,2)*(0.-tr(i,k,it))
	658	!
	659	! if(it.eq.3) write(*,'(I2,1X,a,e20.12,84X,a,e20.12)')k,' 1 det',&
	660	! (-0.+mp(i,2))(qDi(i,k,it)-tr(i,k,it))pdtimeRG/dxpres(i,k)+&
	661	! mp(i,2)(0.-tr(i,k,it))pdtimeRG/dxpres(i,k),&
	662	! 'mp',mp(i,k)
	663	ELSE ! mp(2) = 0 = mp(1) pas de descente insaturee, rien ne se passe s'il ne pleut pas, sinon pluie->env
	664	if(NO_precip(i,1)) THEN
	665	uscavtrac=0.
	666	ELSE
	667	uscavtrac=-imp(i,k)tr(i,k,it)dxpres(i,k)/RG+evap(i,k)qPr(i,k,it)dxpres(i,k)/RG
	668	ENDIF
	669	! if(it.eq.3) write(,'(I2,1X,a,2X,e20.12,82X,a,e20.12)')k,'1 P env incloud',uscavtracpdtimeRG/dxpres(i,k), 'mp',mp(i,k)
	670	ENDIF
	671
	672	ELSE ! k > INB au-dessus du nuage
	673	uscavtrac=0.
	674	ENDIF
	675
	676	! ===== tendances finales ======
	677	trsptd(i,k,it)=trsptrac*pdtimeRG/dxpres(i,k) ! td transport sans eau dans courants satures
	678	dtrSscav(i,k,it)=scavtrac*pdtimeRG/dxpres(i,k) ! td du lessivage dans courants satures
	679	dtrUscav(i,k,it)=uscavtrac*pdtimeRG/dxpres(i,k) ! td courant insat
	680	dtrsat(i,k,it)=(trsptrac+scavtrac)*pdtimeRG/dxpres(i,k) ! td courant sat
	681	dtrcv(i,k,it)=(trsptrac+scavtrac+uscavtrac)*pdtimeRG/dxpres(i,k)!dtrsat(i,k,it)+dtrUscav(i,k,it) td conv
	682	!!!!!!
	683	dtrcvMA(i,k,it)=tdcvMA*pdtimeRG/dxpres(i,k) ! MA tendance convection
	684	ENDDO
	685	ENDDO
	686
	687	! test de conservation du traceur
	688	!print*,"_____________________________________________________________"
	689	!print*," "
	690	! conserv=0.
	691	! conservMA=0.
	692	! DO k= klev-1,1,-1
	693	! DO i=1, klon
	694	! conserv=conserv+dtrcv(i,k,it)* &
	695	! (paprs(i,k)-paprs(i,k+1))/RG
	696	! conservMA=conservMA+dtrcvMA(i,k,it)* &
	697	! (paprs(i,k)-paprs(i,k+1))/RG
	698	!
	699	! if(it.eq.3) write(*,'(I2,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12)') k,&
	700	! 'MA td ',dtrcvMA(i,k,it)*dxpres(i,k)/RG,&
	701	! ' td',dtrcv(i,k,it)*dxpres(i,k)/RG,' conservMA ',conservMA,'conserv ',conserv
	702	!!
	703	! ENDDO
	704	! ENDDO
	705	! if(it.eq.3) print *,'it',it,'conserv ',conserv,'conservMA ',conservMA
	706
	707	END SUBROUTINE cvltr_scav

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

Download in other formats:

Original Format