Context Navigation

cvltr_scav.f90 @ 5274

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

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

Download in other formats:

Original Format