Context Navigation

cvltr_spl.F90 @ 4999

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

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

Download in other formats:

Original Format