Context Navigation

cvltr_scav.F90 @ 3814

Last change on this file since 3814 was 3814, checked in by ymipsl, 10 years ago

remove all dynamic dependency in LMDZ physics except for the include "dimensions.h"

File size: 28.1 KB

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

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

Download in other formats:

Original Format