Context Navigation

cvltr_scav.f90 @ 5450

Last change on this file since 5450 was 5447, checked in by jyg, 5 days ago
output the convective wet deposit of tracers
File size: 29.1 KB

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

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

Download in other formats:

Original Format