Context Navigation

cvltr.F90 @ 2136

Last change on this file since 2136 was 2007, checked in by fhourdin, 11 years ago

Modification pour garantir la conservation des espèces traces.

Modification for tracer conservation

Jean-Yves Grandpeix

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)
	12	USE IOIPSL
[1191]	13	USE dimphy
[1742]	14	USE infotrac, ONLY : nbtr,tname
[1191]	15	IMPLICIT NONE
	16	!=====================================================================
	17	! Objet : convection des traceurs / KE
	18	! Auteurs: M-A Filiberti and J-Y Grandpeix
[1742]	19	! modifiee par R Pilon : lessivage des traceurs / KE
[1191]	20	!=====================================================================
[619]	21
[1191]	22	include "YOMCST.h"
[1742]	23	include "YOECUMF.h"
	24	include "conema3.h"
[1191]	25
	26	! Entree
	27	REAL,INTENT(IN) :: pdtime
	28	REAL,DIMENSION(klon,klev),INTENT(IN) :: da
	29	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi
[1742]	30	! RomP
	31	REAL,DIMENSION(klon,klev),INTENT(IN) :: d1a,dam ! matrices pour simplifier
	32	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi2 ! l'ecriture des tendances
	33	!
	34	REAL,DIMENSION(klon,klev),INTENT(IN) :: mpIN
	35	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression aux 1/2 couches (bas en haut)
	36	! REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression aux 1/2 couches (bas en haut)
	37	REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr ! q de traceur (bas en haut)
	38	INTEGER,INTENT(IN) :: it
[1191]	39	REAL,DIMENSION(klon,klev),INTENT(IN) :: upd ! saturated updraft mass flux
	40	REAL,DIMENSION(klon,klev),INTENT(IN) :: dnd ! saturated downdraft mass flux
[1742]	41	!
	42	REAL,DIMENSION(klon,klev),INTENT(IN) :: wdtrainA ! masses precipitantes de l'asc adiab
	43	REAL,DIMENSION(klon,klev),INTENT(IN) :: wdtrainM ! masses precipitantes des melanges
	44	REAL,DIMENSION(klon,klev),INTENT(IN) :: pmflxrIN ! vprecip: eau
	45	REAL,DIMENSION(klon,klev),INTENT(IN) :: pmflxsIN ! vprecip: neige
	46	REAL,DIMENSION(klon,klev),INTENT(IN) :: ev ! evaporation cv30_routine
	47	REAL,DIMENSION(klon,klev),INTENT(IN) :: epIN
	48	REAL,DIMENSION(klon,klev),INTENT(IN) :: te
	49	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: sij ! fraction dair de lenv
[2007]	50	REAL,DIMENSION(klon,klev),INTENT(IN) :: wght_cvfd ! weights of the layers feeding convection
[1742]	51	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: elij ! contenu en eau condensée spécifique/conc deau condensée massique
	52	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: epmlmMm ! eau condensee precipitee dans mel masse dair sat
	53	REAL,DIMENSION(klon,klev),INTENT(IN) :: eplaMm ! eau condensee precipitee dans aa masse dair sat
[1191]	54
[1742]	55	REAL,DIMENSION(klon,klev),INTENT(IN) :: clw ! contenu en eau condensée dans lasc adiab
	56	REAL,DIMENSION(klon),INTENT(IN) :: sigd
	57	INTEGER,DIMENSION(klon),INTENT(IN) :: icb,inb
[1191]	58	! Sortie
[1742]	59	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrcv ! tendance totale (bas en haut)
	60	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrcvMA ! M-A Filiberti
	61	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: trsptd ! tendance du transport
	62	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrSscav ! tendance du lessivage courant sat
	63	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrsat ! tendance trsp+sat scav
	64	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrUscav ! tendance du lessivage courant unsat
	65	!
	66	! Variables locales
	67	INTEGER :: i,j,k
	68	REAL,DIMENSION(klon,klev) :: dxpres ! difference de pression entre niveau (j+1) et (j)
	69	REAL :: pdtimeRG ! pas de temps * gravite
	70	! variables pour les courants satures
	71	REAL,DIMENSION(klon,klev,klev) :: zmd
	72	REAL,DIMENSION(klon,klev,klev) :: za
	73	REAL,DIMENSION(klon,klev,nbtr) :: zmfd,zmfa
	74	REAL,DIMENSION(klon,klev,nbtr) :: zmfp,zmfu
[2007]	75	REAL,DIMENSION(klon,nbtr) :: qfeed ! tracer concentration feeding convection
[1191]	76
[1742]	77	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfd1a
	78	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfdam
	79	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfphi2
[1191]	80
[1742]	81	! RomP ! les variables sont nettoyees des valeurs aberrantes
	82	REAL,DIMENSION(klon,klev) :: Pa, Pm ! pluie AA et mélanges, var temporaire
	83	REAL,DIMENSION(klon,klev) :: pmflxs,pmflxr ! pmflxrIN,pmflxsIN sans valeur aberante
	84	REAL,DIMENSION(klon,klev) :: mp ! flux de masse
	85	REAL,DIMENSION(klon,klev) :: ep ! fraction d'eau convertie en precipitation
	86	REAL,DIMENSION(klon,klev) :: evap ! evaporation : variable temporaire
	87	REAL,DIMENSION(klon,klev) :: rho !environmental density
	88
	89	REAL,DIMENSION(klon,klev) :: kappa ! denominateur du au calcul de la matrice
	90	! pour obtenir qd et qp
	91	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qTrdi ! traceurs descente air insature transport MA
	92	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qDi ! traceurs descente insaturees
	93	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qPr ! traceurs colonne precipitante
	94	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qPa ! traceurs dans les precip issues lasc. adiab.
	95	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qMel ! traceurs dans les precip issues des melanges
	96	REAL,DIMENSION(klon,klev,nbtr) :: qMeltmp ! variable temporaire
	97	REAL,DIMENSION(klon,klev,nbtr) :: qpmMint
	98	REAL,DIMENSION(klon,klev),INTENT(OUT) :: Mint
	99	! tendances
	100	REAL :: tdcvMA ! terme de transport de traceur (schema Marie Angele)
	101	REAL :: trsptrac ! terme de transport de traceur par l'air
	102	REAL :: scavtrac ! terme de lessivage courant sature
	103	REAL :: uscavtrac ! terme de lessivage courant insature
	104	! impaction
	105	!!! Correction apres discussion Romain P. / Olivier B.
	106	!!! REAL,PARAMETER :: rdrop=2.5e-3 ! rayon des gouttes d'eau
	107	REAL,PARAMETER :: rdrop=1.e-3 ! rayon des gouttes d'eau
	108	!!!
	109	REAL,DIMENSION(klon,klev) :: imp ! coefficient d'impaction
	110	! parametres lessivage
	111	REAL :: ccntrAA_coef ! \alpha_a : fract aerosols de l'AA convertis en CCN
	112	REAL :: ccntrENV_coef ! \beta_m : fract aerosols de l'env convertis en CCN
	113	REAL :: coefcoli ! coefficient de collision des gouttes sur les aerosols
	114	!
	115	LOGICAL,DIMENSION(klon,klev) :: NO_precip
	116	! LOGICAL :: scavON
	117	! var tmp tests
	118	REAL :: conserv
	119	real :: conservMA
	120
	121	! coefficient lessivage
	122	ccntrAA_coef = 0.
	123	ccntrENV_coef = 0.
	124	coefcoli = 0.
	125
[1759]	126	!$OMP MASTER
[1742]	127	call getin('ccntrAA_coef',ccntrAA_coef)
	128	call getin('ccntrENV_coef',ccntrENV_coef)
	129	call getin('coefcoli',coefcoli)
[1759]	130	!$OMP END MASTER
	131	!$OMP BARRIER
[1742]	132	print*,'cvltr coef lessivage convectif', ccntrAA_coef,ccntrENV_coef,coefcoli
	133
	134	! scavON=.TRUE.
	135	! if(scavON) then
	136	! ccntrAA_coef = 1.
	137	! ccntrENV_coef = 1.
	138	! coefcoli = 1.
	139	! else
	140	! ccntrAA_coef = 0.
	141	! ccntrENV_coef = 0.
	142	! coefcoli = 0.
	143	! endif
	144
	145	! ======================================================
	146	! calcul de l'impaction
	147	! ======================================================
	148	!initialisation
	149	do j=1,klev
	150	do i=1,klon
	151	imp(i,j)=0.
	152	enddo
	153	enddo
	154	! impaction sur la surface de la colonne de la descente insaturee
	155	! On prend la moyenne des precip entre le niveau i+1 et i
	156	! I=3/4* (P(1+1)+P(i))/2 / (sigdrrho_l)
	157	! 1000kg/m3= densité de l'eau
	158	! 0.75e-3 = 3/4 /1000
	159	! Par la suite, I est tout le temps multiplié par sig_d pour avoir l'impaction sur la surface de la maille
	160	! on le néglige ici pour simplifier le code
	161	do j=1,klev-1
	162	do i=1,klon
	163	imp(i,j) = coefcoli0.75e-3/rdrop &
	164	0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j))
	165	! rho(i,j)=pplay(i,j)/(rd*te(i,j))
	166	enddo
	167	enddo
	168	!
	169	! initialisation pour flux de traceurs, td et autre
	170	trsptrac = 0.
	171	scavtrac = 0.
	172	uscavtrac = 0.
[2007]	173	qfeed(:,it) = 0. !RL
[1250]	174	DO j=1,klev
[1742]	175	DO i=1,klon
	176	zmfd(i,j,it)=0.
	177	zmfa(i,j,it)=0.
	178	zmfu(i,j,it)=0.
	179	zmfp(i,j,it)=0.
	180	zmfphi2(i,j,it)=0.
	181	zmfd1a(i,j,it)=0.
	182	zmfdam(i,j,it)=0.
	183	qDi(i,j,it)=0.
	184	qPr(i,j,it)=0.
	185	qPa(i,j,it)=0.
	186	qMel(i,j,it)=0.
	187	qMeltmp(i,j,it)=0.
	188	qTrdi(i,j,it)=0.
	189	kappa(i,j)=0.
	190	trsptd(i,j,it)=0.
	191	dtrsat(i,j,it)=0.
	192	dtrSscav(i,j,it)=0.
	193	dtrUscav(i,j,it)=0.
	194	dtrcv(i,j,it)=0.
	195	dtrcvMA(i,j,it)=0.
	196	evap(i,j)=0.
	197	dxpres(i,j)=0.
	198	qpmMint(i,j,it)=0.
	199	Mint(i,j)=0.
	200	END DO
[1250]	201	END DO
[1742]	202
	203	! suppression des valeurs très faibles (~1e-320)
	204	! multiplication de levaporation pour lavoir par unite de temps
	205	! et par unite de surface de la maille
	206	! -> cv30_unsat : evap : masse evaporee/s/(m2 de la descente)
	207	DO j=1,klev
	208	DO i=1,klon
	209	if(ev(i,j).lt.1.e-16) then
	210	evap(i,j)=0.
	211	else
	212	evap(i,j)=ev(i,j)*sigd(i)
	213	endif
	214	END DO
[1250]	215	END DO
[1742]	216
[1250]	217	DO j=1,klev
[1742]	218	DO i=1,klon
	219	if(j.lt.klev) then
	220	if(epIN(i,j).lt.1.e-32) then
	221	ep(i,j)=0.
	222	else
	223	ep(i,j)=epIN(i,j)
	224	endif
	225	else
	226	ep(i,j)=epmax
	227	endif
	228	if(mpIN(i,j).lt.1.e-32) then
	229	mp(i,j)=0.
	230	else
	231	mp(i,j)=mpIN(i,j)
	232	endif
	233	if(pmflxsIN(i,j).lt.1.e-32) then
	234	pmflxs(i,j)=0.
	235	else
	236	pmflxs(i,j)=pmflxsIN(i,j)
	237	endif
	238	if(pmflxrIN(i,j).lt.1.e-32) then
	239	pmflxr(i,j)=0.
	240	else
	241	pmflxr(i,j)=pmflxrIN(i,j)
	242	endif
	243	if(wdtrainA(i,j).lt.1.e-32) then
	244	Pa(i,j)=0.
	245	else
	246	Pa(i,j)=wdtrainA(i,j)
	247	endif
	248	if(wdtrainM(i,j).lt.1.e-32) then
	249	Pm(i,j)=0.
	250	else
	251	Pm(i,j)=wdtrainM(i,j)
	252	endif
	253	END DO
[1250]	254	END DO
[1742]	255
	256	!==========================================
	257	DO j = klev-1,1,-1
	258	DO i = 1,klon
	259	NO_precip(i,j) = (pmflxr(i,j+1)+pmflxs(i,j+1)).lt.1.e-10&
	260	.and.Pa(i,j).lt.1.e-10.and.Pm(i,j).lt.1.e-10
	261	END DO
[1250]	262	END DO
[1742]	263
	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
[1250]	275	END DO
[1191]	276	! calcul de la matrice d echange
	277	! matrice de distribution de la masse entrainee en k
[1742]	278	! commmentaire RomP : mp > 0
[1250]	279	DO k=1,klev-1
[1191]	280	DO i=1,klon
[1742]	281	zmd(i,k,k)=max(0.,mp(i,k)-mp(i,k+1)) ! ~ mk(k)
[1191]	282	END DO
	283	END DO
	284	DO k=2,klev
	285	DO j=k-1,1,-1
	286	DO i=1,klon
[1742]	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)
[1191]	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
[1742]	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	!!!!!
[1191]	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	!
[1742]	322	! tendance courants insatures ! sans lessivage ancien schema
	323	!
[1191]	324	DO k=1,klev
	325	DO j=1,klev
	326	DO i=1,klon
[1742]	327	zmfd(i,j,it)=zmfd(i,j,it)+za(i,j,k)*(tr(i,k,it)-tr(i,j,it))
[1191]	328	END DO
	329	END DO
	330	END DO
	331	!
	332	! =========================================
[1742]	333	! calcul des tendances liees aux courants satures j <-> z ; k <-> z'
[1191]	334	! =========================================
[2007]	335	!
	336	!RL
	337	! Feeding concentrations
[1191]	338	DO j=1,klev
	339	DO i=1,klon
[2007]	340	qfeed(i,it)=qfeed(i,it)+wght_cvfd(i,j)*tr(i,j,it)
[1191]	341	END DO
	342	END DO
[2007]	343	!RL
	344	!
	345	DO j=1,klev
	346	DO i=1,klon
	347	!RL
	348	!! zmfa(i,j,it)=da(i,j)*(tr(i,1,it)-tr(i,j,it)) ! da
	349	zmfa(i,j,it)=da(i,j)*(qfeed(i,it)-tr(i,j,it)) ! da
	350	!RL
	351	END DO
	352	END DO
	353	!
[1191]	354	DO k=1,klev
	355	DO j=1,klev
	356	DO i=1,klon
[1742]	357	zmfp(i,j,it)=zmfp(i,j,it)+phi(i,j,k)*(tr(i,k,it)-tr(i,j,it)) ! phi
[1191]	358	END DO
	359	END DO
	360	END DO
[1742]	361	! RomP ajout des matrices liees au lessivage
	362	DO j=1,klev
	363	DO i=1,klon
	364	zmfd1a(i,j,it)=d1a(i,j)*tr(i,1,it) ! da1
	365	zmfdam(i,j,it)=dam(i,j)*tr(i,1,it) ! dam
	366	END DO
	367	END DO
	368	DO k=1,klev
	369	DO j=1,klev
	370	DO i=1,klon
	371	zmfphi2(i,j,it)=zmfphi2(i,j,it)+phi2(i,j,k)*tr(i,k,it) ! psi
	372	END DO
	373	END DO
	374	END DO
[1191]	375	DO j=1,klev-1
	376	DO i=1,klon
[1742]	377	zmfu(i,j,it)=max(0.,upd(i,j+1)+dnd(i,j+1))*(tr(i,j+1,it)-tr(i,j,it))
[1191]	378	END DO
	379	END DO
	380	DO j=2,klev
	381	DO i=1,klon
[1742]	382	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]	383	END DO
	384	END DO
[1742]	385	! ===================================================
	386	! calcul des tendances liees aux courants insatures
	387	! ===================================================
	388	! pression
[1191]	389	DO k=1, klev
	390	DO i=1, klon
[1742]	391	dxpres(i,k)=paprs(i,k)-paprs(i,k+1)
[1250]	392	ENDDO
	393	ENDDO
	394	pdtimeRG=pdtime*RG
[1742]	395
	396	! q_pa et q_pm traceurs issues des courants satures se retrouvant dans les precipitations
	397	DO j=1,klev
	398	DO i=1,klon
	399	if(j.ge.icb(i).and.j.le.inb(i)) then
	400	if(clw(i,j).gt.1.e-16) then
	401	qPa(i,j,it)=ccntrAA_coef*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	qMeltmp(i,j,it)=((1-ep(i,k))ccntrAA_coeftr(i,1,it)&
	419	*(1.-sij(i,k,j)) +ccntrENV_coef&
	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	scavtrac=-ccntrAA_coef*zmfd1a(i,k,it)&
	613	-zmfphi2(i,k,it)*ccntrENV_coef&
	614	-zmfdam(i,k,it)*ccntrAA_coef
	615	! lessivage courants insatures
	616	if(k.le.inb(i).and.k.gt.1) then ! tendances dans le nuage
	617	!------------------------------------------------------------- detrainement
	618	if(mp(i,k+1).gt.mp(i,k).and.mp(i,k+1).gt.1.e-10) then
	619	uscavtrac= (-mp(i,k)+mp(i,k+1))*(qDi(i,k,it)-tr(i,k,it))&
	620	+ mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))
	621	!
	622	! if(it.eq.3) write(*,'(I2,1X,a,5X,e20.12,82X,a,e20.12)')k,' det incloud',&
	623	! (-mp(i,k)+mp(i,k+1))(qDi(i,k,it)-tr(i,k,it))pdtimeRG/dxpres(i,k)+&
	624	! mp(i,k)(tr(i,k-1,it)-tr(i,k,it))pdtimeRG/dxpres(i,k),&
	625	! 'mp',mp(i,k)
	626	!------------------------------------------------------------- entrainement
	627	elseif(mp(i,k).gt.mp(i,k+1).and.mp(i,k).gt.1.e-10) then
	628	uscavtrac= mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))
	629	!
	630	! 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)
	631	!=!------------------------------------------------------------- end ent/det
	632	else ! mp(i,k+1)=0. et mp(i,k)=0. pluie directement sur l environnement
	633
	634	if(NO_precip(i,k)) then
	635	uscavtrac=0.
	636	! 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)
	637	else
	638	uscavtrac=-imp(i,k)tr(i,k,it)dxpres(i,k)/RG+evap(i,k)qPr(i,k,it)dxpres(i,k)/RG
	639	! 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)
	640	endif
	641	endif ! mp/det/ent
	642	!------------------------------------------------------------- premiere couche
	643	elseif(k.eq.1) then ! mp(1)=0.
	644	if(mp(i,2).gt.1.e-10) then !detrainement
	645	uscavtrac= (-0.+mp(i,2))*(qDi(i,k,it)-tr(i,k,it)) !&
	646	! + mp(i,2)*(0.-tr(i,k,it))
	647	!
	648	! if(it.eq.3) write(*,'(I2,1X,a,e20.12,84X,a,e20.12)')k,' 1 det',&
	649	! (-0.+mp(i,2))(qDi(i,k,it)-tr(i,k,it))pdtimeRG/dxpres(i,k)+&
	650	! mp(i,2)(0.-tr(i,k,it))pdtimeRG/dxpres(i,k),&
	651	! 'mp',mp(i,k)
	652	else ! mp(2) = 0 = mp(1) pas de descente insaturee, rien ne se passe s'il ne pleut pas, sinon pluie->env
	653	if(NO_precip(i,1)) then
	654	uscavtrac=0.
	655	else
	656	uscavtrac=-imp(i,k)tr(i,k,it)dxpres(i,k)/RG+evap(i,k)qPr(i,k,it)dxpres(i,k)/RG
	657	endif
	658	! 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)
	659	endif
	660
	661	else ! k > INB au-dessus du nuage
	662	uscavtrac=0.
	663	endif
	664
	665	! ===== tendances finales ======
	666	trsptd(i,k,it)=trsptrac*pdtimeRG/dxpres(i,k) ! td transport sans eau dans courants satures
	667	dtrSscav(i,k,it)=scavtrac*pdtimeRG/dxpres(i,k) ! td du lessivage dans courants satures
	668	dtrUscav(i,k,it)=uscavtrac*pdtimeRG/dxpres(i,k) ! td courant insat
	669	dtrsat(i,k,it)=(trsptrac+scavtrac)*pdtimeRG/dxpres(i,k) ! td courant sat
	670	dtrcv(i,k,it)=(trsptrac+scavtrac+uscavtrac)*pdtimeRG/dxpres(i,k)!dtrsat(i,k,it)+dtrUscav(i,k,it) td conv
	671	!!!!!!
	672	dtrcvMA(i,k,it)=tdcvMA*pdtimeRG/dxpres(i,k) ! MA tendance convection
[1191]	673	ENDDO
	674	ENDDO
	675
	676	! test de conservation du traceur
[1742]	677	!print*,"_____________________________________________________________"
	678	!print*," "
[1191]	679	! conserv=0.
[1742]	680	! conservMA=0.
	681	! DO k= klev-1,1,-1
[1191]	682	! DO i=1, klon
[1742]	683	! conserv=conserv+dtrcv(i,k,it)* &
[1191]	684	! (paprs(i,k)-paprs(i,k+1))/RG
[1742]	685	! conservMA=conservMA+dtrcvMA(i,k,it)* &
	686	! (paprs(i,k)-paprs(i,k+1))/RG
	687	!
	688	! if(it.eq.3) write(*,'(I2,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12)') k,&
	689	! 'MA td ',dtrcvMA(i,k,it)*dxpres(i,k)/RG,&
	690	! ' td',dtrcv(i,k,it)*dxpres(i,k)/RG,' conservMA ',conservMA,'conserv ',conserv
	691	!!
[1191]	692	! ENDDO
	693	! ENDDO
[1742]	694	! if(it.eq.3) print *,'it',it,'conserv ',conserv,'conservMA ',conservMA
	695
[1191]	696	END SUBROUTINE cvltr

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

Download in other formats:

Original Format