Context Navigation

cvltr_spl.f90 @ 5300

Last change on this file since 5300 was 5292, checked in by abarral, 5 weeks ago
Move academic.h chem.h chem_spla.h to module
File size: 30.2 KB

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

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

Download in other formats:

Original Format