Context Navigation

cvltr_scav.f90 @ 5283

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

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

Download in other formats:

Original Format