Context Navigation

cvltr_scav.F90 @ 3230

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

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

Download in other formats:

Original Format