Context Navigation

cvltr.f90 @ 5424

Last change on this file since 5424 was 5289, checked in by abarral, 7 weeks ago
Turn YOECUMF.h into a module Fix USE in fxy_new_mod_h.f90
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 25.6 KB

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

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

Download in other formats:

Original Format