Context Navigation

cvltr.F90 @ 1934

Last change on this file since 1934 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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.1 KB

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

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

Download in other formats:

Original Format