Context Navigation

cvltr_scav.f90 @ 5559

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

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

Download in other formats:

Original Format