Context Navigation

cvltr_scav.f90 @ 5447

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

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

Download in other formats:

Original Format