Context Navigation

cvltr_scav.f90 @ 5441

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

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

Download in other formats:

Original Format