Context Navigation

cvltr.F90 @ 1910

Last change on this file since 1910 was 1910, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes r1860:1909 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 25.1 KB

Line
1	!
2	! $Id $
3	!
4	SUBROUTINE cvltr(pdtime, da, phi,phi2,d1a,dam, mpIN,epIN, &
5	sigd,sij,clw,elij,epmlmMm,eplaMm, &
6	pmflxrIN,pmflxsIN,ev,te,wdtrainA,wdtrainM, &
7	paprs,it,tr,upd,dnd,inb,icb, &
8	dtrcv,trsptd,dtrSscav,dtrsat,dtrUscav,qDi,qPr, &
9	qPa,qMel,qTrdi,dtrcvMA,Mint, &
10	zmfd1a,zmfphi2,zmfdam)
11	USE IOIPSL
12	USE dimphy
13	USE infotrac, ONLY : nbtr,tname
14	IMPLICIT NONE
15	!=====================================================================
16	! Objet : convection des traceurs / KE
17	! Auteurs: M-A Filiberti and J-Y Grandpeix
18	! modifiee par R Pilon : lessivage des traceurs / KE
19	!=====================================================================
20
21	include "YOMCST.h"
22	include "YOECUMF.h"
23	include "conema3.h"
24
25	! Entree
26	REAL,INTENT(IN) :: pdtime
27	REAL,DIMENSION(klon,klev),INTENT(IN) :: da
28	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi
29	! RomP
30	REAL,DIMENSION(klon,klev),INTENT(IN) :: d1a,dam ! matrices pour simplifier
31	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi2 ! l'ecriture des tendances
32	!
33	REAL,DIMENSION(klon,klev),INTENT(IN) :: mpIN
34	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression aux 1/2 couches (bas en haut)
35	! REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression aux 1/2 couches (bas en haut)
36	REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr ! q de traceur (bas en haut)
37	INTEGER,INTENT(IN) :: it
38	REAL,DIMENSION(klon,klev),INTENT(IN) :: upd ! saturated updraft mass flux
39	REAL,DIMENSION(klon,klev),INTENT(IN) :: dnd ! saturated downdraft mass flux
40	!
41	REAL,DIMENSION(klon,klev),INTENT(IN) :: wdtrainA ! masses precipitantes de l'asc adiab
42	REAL,DIMENSION(klon,klev),INTENT(IN) :: wdtrainM ! masses precipitantes des melanges
43	REAL,DIMENSION(klon,klev),INTENT(IN) :: pmflxrIN ! vprecip: eau
44	REAL,DIMENSION(klon,klev),INTENT(IN) :: pmflxsIN ! vprecip: neige
45	REAL,DIMENSION(klon,klev),INTENT(IN) :: ev ! evaporation cv30_routine
46	REAL,DIMENSION(klon,klev),INTENT(IN) :: epIN
47	REAL,DIMENSION(klon,klev),INTENT(IN) :: te
48	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: sij ! fraction dair de lenv
49	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: elij ! contenu en eau condensée spécifique/conc deau condensée massique
50	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: epmlmMm ! eau condensee precipitee dans mel masse dair sat
51	REAL,DIMENSION(klon,klev),INTENT(IN) :: eplaMm ! eau condensee precipitee dans aa masse dair sat
52
53	REAL,DIMENSION(klon,klev),INTENT(IN) :: clw ! contenu en eau condensée dans lasc adiab
54	REAL,DIMENSION(klon),INTENT(IN) :: sigd
55	INTEGER,DIMENSION(klon),INTENT(IN) :: icb,inb
56	! Sortie
57	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrcv ! tendance totale (bas en haut)
58	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrcvMA ! M-A Filiberti
59	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: trsptd ! tendance du transport
60	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrSscav ! tendance du lessivage courant sat
61	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrsat ! tendance trsp+sat scav
62	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrUscav ! tendance du lessivage courant unsat
63	!
64	! Variables locales
65	INTEGER :: i,j,k
66	REAL,DIMENSION(klon,klev) :: dxpres ! difference de pression entre niveau (j+1) et (j)
67	REAL :: pdtimeRG ! pas de temps * gravite
68	! variables pour les courants satures
69	REAL,DIMENSION(klon,klev,klev) :: zmd
70	REAL,DIMENSION(klon,klev,klev) :: za
71	REAL,DIMENSION(klon,klev,nbtr) :: zmfd,zmfa
72	REAL,DIMENSION(klon,klev,nbtr) :: zmfp,zmfu
73
74	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfd1a
75	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfdam
76	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfphi2
77
78	! RomP ! les variables sont nettoyees des valeurs aberrantes
79	REAL,DIMENSION(klon,klev) :: Pa, Pm ! pluie AA et mélanges, var temporaire
80	REAL,DIMENSION(klon,klev) :: pmflxs,pmflxr ! pmflxrIN,pmflxsIN sans valeur aberante
81	REAL,DIMENSION(klon,klev) :: mp ! flux de masse
82	REAL,DIMENSION(klon,klev) :: ep ! fraction d'eau convertie en precipitation
83	REAL,DIMENSION(klon,klev) :: evap ! evaporation : variable temporaire
84	REAL,DIMENSION(klon,klev) :: rho !environmental density
85
86	REAL,DIMENSION(klon,klev) :: kappa ! denominateur du au calcul de la matrice
87	! pour obtenir qd et qp
88	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qTrdi ! traceurs descente air insature transport MA
89	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qDi ! traceurs descente insaturees
90	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qPr ! traceurs colonne precipitante
91	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qPa ! traceurs dans les precip issues lasc. adiab.
92	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qMel ! traceurs dans les precip issues des melanges
93	REAL,DIMENSION(klon,klev,nbtr) :: qMeltmp ! variable temporaire
94	REAL,DIMENSION(klon,klev,nbtr) :: qpmMint
95	REAL,DIMENSION(klon,klev),INTENT(OUT) :: Mint
96	! tendances
97	REAL :: tdcvMA ! terme de transport de traceur (schema Marie Angele)
98	REAL :: trsptrac ! terme de transport de traceur par l'air
99	REAL :: scavtrac ! terme de lessivage courant sature
100	REAL :: uscavtrac ! terme de lessivage courant insature
101	! impaction
102	!!! Correction apres discussion Romain P. / Olivier B.
103	!!! REAL,PARAMETER :: rdrop=2.5e-3 ! rayon des gouttes d'eau
104	REAL,PARAMETER :: rdrop=1.e-3 ! rayon des gouttes d'eau
105	!!!
106	REAL,DIMENSION(klon,klev) :: imp ! coefficient d'impaction
107	! parametres lessivage
108	REAL :: ccntrAA_coef ! \alpha_a : fract aerosols de l'AA convertis en CCN
109	REAL :: ccntrENV_coef ! \beta_m : fract aerosols de l'env convertis en CCN
110	REAL :: coefcoli ! coefficient de collision des gouttes sur les aerosols
111	!
112	LOGICAL,DIMENSION(klon,klev) :: NO_precip
113	! LOGICAL :: scavON
114	! var tmp tests
115	REAL :: conserv
116	real :: conservMA
117
118	! coefficient lessivage
119	ccntrAA_coef = 0.
120	ccntrENV_coef = 0.
121	coefcoli = 0.
122
123	!$OMP MASTER
124	call getin('ccntrAA_coef',ccntrAA_coef)
125	call getin('ccntrENV_coef',ccntrENV_coef)
126	call getin('coefcoli',coefcoli)
127	!$OMP END MASTER
128	!$OMP BARRIER
129	print*,'cvltr coef lessivage convectif', ccntrAA_coef,ccntrENV_coef,coefcoli
130
131	! scavON=.TRUE.
132	! if(scavON) then
133	! ccntrAA_coef = 1.
134	! ccntrENV_coef = 1.
135	! coefcoli = 1.
136	! else
137	! ccntrAA_coef = 0.
138	! ccntrENV_coef = 0.
139	! coefcoli = 0.
140	! endif
141
142	! ======================================================
143	! calcul de l'impaction
144	! ======================================================
145	!initialisation
146	do j=1,klev
147	do i=1,klon
148	imp(i,j)=0.
149	enddo
150	enddo
151	! impaction sur la surface de la colonne de la descente insaturee
152	! On prend la moyenne des precip entre le niveau i+1 et i
153	! I=3/4* (P(1+1)+P(i))/2 / (sigdrrho_l)
154	! 1000kg/m3= densité de l'eau
155	! 0.75e-3 = 3/4 /1000
156	! Par la suite, I est tout le temps multiplié par sig_d pour avoir l'impaction sur la surface de la maille
157	! on le néglige ici pour simplifier le code
158	do j=1,klev-1
159	do i=1,klon
160	imp(i,j) = coefcoli0.75e-3/rdrop &
161	0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j))
162	! rho(i,j)=pplay(i,j)/(rd*te(i,j))
163	enddo
164	enddo
165	!
166	! initialisation pour flux de traceurs, td et autre
167	trsptrac = 0.
168	scavtrac = 0.
169	uscavtrac = 0.
170
171	DO j=1,klev
172	DO i=1,klon
173	zmfd(i,j,it)=0.
174	zmfa(i,j,it)=0.
175	zmfu(i,j,it)=0.
176	zmfp(i,j,it)=0.
177	zmfphi2(i,j,it)=0.
178	zmfd1a(i,j,it)=0.
179	zmfdam(i,j,it)=0.
180	qDi(i,j,it)=0.
181	qPr(i,j,it)=0.
182	qPa(i,j,it)=0.
183	qMel(i,j,it)=0.
184	qMeltmp(i,j,it)=0.
185	qTrdi(i,j,it)=0.
186	kappa(i,j)=0.
187	trsptd(i,j,it)=0.
188	dtrsat(i,j,it)=0.
189	dtrSscav(i,j,it)=0.
190	dtrUscav(i,j,it)=0.
191	dtrcv(i,j,it)=0.
192	dtrcvMA(i,j,it)=0.
193	evap(i,j)=0.
194	dxpres(i,j)=0.
195	qpmMint(i,j,it)=0.
196	Mint(i,j)=0.
197	END DO
198	END DO
199
200	! suppression des valeurs très faibles (~1e-320)
201	! multiplication de levaporation pour lavoir par unite de temps
202	! et par unite de surface de la maille
203	! -> cv30_unsat : evap : masse evaporee/s/(m2 de la descente)
204	DO j=1,klev
205	DO i=1,klon
206	if(ev(i,j).lt.1.e-16) then
207	evap(i,j)=0.
208	else
209	evap(i,j)=ev(i,j)*sigd(i)
210	endif
211	END DO
212	END DO
213
214	DO j=1,klev
215	DO i=1,klon
216	if(j.lt.klev) then
217	if(epIN(i,j).lt.1.e-32) then
218	ep(i,j)=0.
219	else
220	ep(i,j)=epIN(i,j)
221	endif
222	else
223	ep(i,j)=epmax
224	endif
225	if(mpIN(i,j).lt.1.e-32) then
226	mp(i,j)=0.
227	else
228	mp(i,j)=mpIN(i,j)
229	endif
230	if(pmflxsIN(i,j).lt.1.e-32) then
231	pmflxs(i,j)=0.
232	else
233	pmflxs(i,j)=pmflxsIN(i,j)
234	endif
235	if(pmflxrIN(i,j).lt.1.e-32) then
236	pmflxr(i,j)=0.
237	else
238	pmflxr(i,j)=pmflxrIN(i,j)
239	endif
240	if(wdtrainA(i,j).lt.1.e-32) then
241	Pa(i,j)=0.
242	else
243	Pa(i,j)=wdtrainA(i,j)
244	endif
245	if(wdtrainM(i,j).lt.1.e-32) then
246	Pm(i,j)=0.
247	else
248	Pm(i,j)=wdtrainM(i,j)
249	endif
250	END DO
251	END DO
252
253	!==========================================
254	DO j = klev-1,1,-1
255	DO i = 1,klon
256	NO_precip(i,j) = (pmflxr(i,j+1)+pmflxs(i,j+1)).lt.1.e-10&
257	.and.Pa(i,j).lt.1.e-10.and.Pm(i,j).lt.1.e-10
258	END DO
259	END DO
260
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	DO j=1,klev
333	DO i=1,klon
334	zmfa(i,j,it)=da(i,j)*(tr(i,1,it)-tr(i,j,it)) ! da
335	END DO
336	END DO
337	DO k=1,klev
338	DO j=1,klev
339	DO i=1,klon
340	zmfp(i,j,it)=zmfp(i,j,it)+phi(i,j,k)*(tr(i,k,it)-tr(i,j,it)) ! phi
341	END DO
342	END DO
343	END DO
344	! RomP ajout des matrices liees au lessivage
345	DO j=1,klev
346	DO i=1,klon
347	zmfd1a(i,j,it)=d1a(i,j)*tr(i,1,it) ! da1
348	zmfdam(i,j,it)=dam(i,j)*tr(i,1,it) ! dam
349	END DO
350	END DO
351	DO k=1,klev
352	DO j=1,klev
353	DO i=1,klon
354	zmfphi2(i,j,it)=zmfphi2(i,j,it)+phi2(i,j,k)*tr(i,k,it) ! psi
355	END DO
356	END DO
357	END DO
358	DO j=1,klev-1
359	DO i=1,klon
360	zmfu(i,j,it)=max(0.,upd(i,j+1)+dnd(i,j+1))*(tr(i,j+1,it)-tr(i,j,it))
361	END DO
362	END DO
363	DO j=2,klev
364	DO i=1,klon
365	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))
366	END DO
367	END DO
368	! ===================================================
369	! calcul des tendances liees aux courants insatures
370	! ===================================================
371	! pression
372	DO k=1, klev
373	DO i=1, klon
374	dxpres(i,k)=paprs(i,k)-paprs(i,k+1)
375	ENDDO
376	ENDDO
377	pdtimeRG=pdtime*RG
378
379	! q_pa et q_pm traceurs issues des courants satures se retrouvant dans les precipitations
380	DO j=1,klev
381	DO i=1,klon
382	if(j.ge.icb(i).and.j.le.inb(i)) then
383	if(clw(i,j).gt.1.e-16) then
384	qPa(i,j,it)=ccntrAA_coef*tr(i,1,it)/clw(i,j)
385	else
386	qPa(i,j,it)=0.
387	endif
388	endif
389	END DO
390	END DO
391
392	! calcul de q_pm en 2 parties :
393	! 1) calcul de sa valeur pour un niveau z' donne
394	! 2) integration sur la verticale sur z'
395	DO j=1,klev
396	DO k=1,j-1
397	DO i=1,klon
398	if(k.ge.icb(i).and.k.le.inb(i).and.&
399	j.le.inb(i)) then
400	if(elij(i,k,j).gt.1.e-16) then
401	qMeltmp(i,j,it)=((1-ep(i,k))ccntrAA_coeftr(i,1,it)&
402	*(1.-sij(i,k,j)) +ccntrENV_coef&
403	tr(i,k,it)sij(i,k,j)) / elij(i,k,j)
404	else
405	qMeltmp(i,j,it)=0.
406	endif
407	qpmMint(i,j,it)=qpmMint(i,j,it) + qMeltmp(i,j,it)*epmlmMm(i,j,k)
408	Mint(i,j)=Mint(i,j) + epmlmMm(i,j,k)
409	endif ! end if dans nuage
410	END DO
411	END DO
412	END DO
413
414	DO j=1,klev
415	DO i=1,klon
416	if(Mint(i,j).gt.1.e-16) then
417	qMel(i,j,it)=qpmMint(i,j,it)/Mint(i,j)
418	else
419	qMel(i,j,it)=0.
420	endif
421	END DO
422	END DO
423
424	! calcul de q_d et q_p traceurs de la descente precipitante
425	DO j=klev-1,1,-1
426	DO i=1,klon
427	if(mp(i,j+1).gt.mp(i,j).and.mp(i,j+1).gt.1.e-10) then ! detrainement
428	kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
429	(-mp(i,j+1)-imp(i,j)/RG*dxpres(i,j))&
430	+ (imp(i,j)/RGdxpres(i,j))(evap(i,j)/RG*dxpres(i,j)))
431
432	elseif(mp(i,j).gt.mp(i,j+1).and.mp(i,j).gt.1.e-10) then! entrainement
433	if(j.eq.1) then
434	kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
435	(-mp(i,2)-imp(i,j)/RG*dxpres(i,j))&
436	+ (imp(i,j)/RGdxpres(i,j))(evap(i,j)/RG*dxpres(i,j)))
437	else
438	kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
439	(-mp(i,j)-imp(i,j)/RG*dxpres(i,j))&
440	+ (imp(i,j)/RGdxpres(i,j))(evap(i,j)/RG*dxpres(i,j)))
441	endif
442	else
443	kappa(i,j)=1.
444	endif
445	ENDDO
446	ENDDO
447
448	DO j=klev-1,1,-1
449	DO i=1,klon
450	if (abs(kappa(i,j)).lt.1.e-25) then !si denominateur nul (il peut y avoir des mp!=0)
451	kappa(i,j)=1.
452	if(j.eq.1) then
453	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
454	elseif(mp(i,j+1).gt.mp(i,j).and.mp(i,j+1).gt.1.e-10) then
455	qDi(i,j,it)=qDi(i,j+1,it)
456	elseif(mp(i,j).gt.mp(i,j+1).and.mp(i,j).gt.1.e-10) then! entrainement
457	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))
458	else ! si mp (i)=0 et mp(j+1)=0
459	qDi(i,j,it)=tr(i,j,it) ! orig 0.
460	endif
461
462	if(NO_precip(i,j)) then
463	qPr(i,j,it)=0.
464	else
465	qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
466	Pa(i,j)qPa(i,j,it)+Pm(i,j)qMel(i,j,it)&
467	+imp(i,j)/RGdxpres(i,j)qDi(i,j,it))/&
468	(pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))
469	endif
470	else ! denominateur non nul
471	kappa(i,j)=1./kappa(i,j)
472	! calcul de qd et qp
473	!!jyg (20130119) correction pour le sommet du nuage
474	!! if(j.ge.inb(i)) then !au-dessus du nuage, sommet inclu
475	if(j.gt.inb(i)) then !au-dessus du nuage
476	qDi(i,j,it)=tr(i,j,it) ! pas de descente => environnement = descente insaturee
477	qPr(i,j,it)=0.
478
479	! vvv premiere couche du modele ou mp(1)=0 ! det tout le temps vvv
480	elseif(j.eq.1) then
481	if(mp(i,2).gt.1.e-10) then !mp(2) non nul -> detrainement (car mp(1) = 0) !ent pas possible
482	if(NO_precip(i,j)) then !pas de precip en (i)
483	qDi(i,j,it)=qDi(i,j+1,it)
484	qPr(i,j,it)=0.
485	else
486	qDi(i,j,it)=kappa(i,j)*(&
487	(-evap(i,j)/RGdxpres(i,j))((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
488	Pa(i,j)qPa(i,j,it)+Pm(i,j)qMel(i,j,it)) +&
489	(pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
490	(-mp(i,j+1)*qDi(i,j+1,it)))
491
492	qPr(i,j,it)=kappa(i,j)*(&
493	(-mp(i,j+1)-imp(i,j)/RGdxpres(i,j))&
494	((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	+(-mp(i,j+1)qDi(i,j+1,it)) (imp(i,j)/RG*dxpres(i,j)))
497	endif
498
499	else !mp(2) nul -> plus de descente insaturee -> pluie agit sur environnement
500	qDi(i,j,it)=tr(i,j,it) ! orig 0.
501	if(NO_precip(i,j)) then
502	qPr(i,j,it)=0.
503	else
504	qPr(i,j,it)=((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	+imp(i,j)/RGdxpres(i,j)tr(i,j,it))/&
507	(pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))
508	endif
509
510	endif !mp(2) nul ou non
511
512	! vvv (j!=1.and.j.lt.inb(i)) en-dessous du sommet nuage vvv
513	else
514	!------------------------------------------------------------- detrainement
515	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
516	if(NO_precip(i,j)) then
517	qDi(i,j,it)=qDi(i,j+1,it)
518	qPr(i,j,it)=0.
519	else
520	qDi(i,j,it)=kappa(i,j)*(&
521	(-evap(i,j)/RGdxpres(i,j))((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	(pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
524	(-mp(i,j+1)*qDi(i,j+1,it)))
525	!
526	qPr(i,j,it)=kappa(i,j)*(&
527	(-mp(i,j+1)-imp(i,j)/RGdxpres(i,j))&
528	((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
529	Pa(i,j)qPa(i,j,it)+Pm(i,j)qMel(i,j,it))&
530	+(-mp(i,j+1)qDi(i,j+1,it)) (imp(i,j)/RG*dxpres(i,j)))
531	endif !precip
532	!------------------------------------------------------------- entrainement
533	elseif(mp(i,j).gt.mp(i,j+1).and.mp(i,j).gt.1.e-10) then
534	if(NO_precip(i,j)) then
535	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))
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)-tr(i,j,it))-mp(i,j)tr(i,j,it)))
543	!
544	qPr(i,j,it)=kappa(i,j)*(&
545	(-mp(i,j)-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)-tr(i,j,it))-mp(i,j)tr(i,j,it))*&
549	(imp(i,j)/RG*dxpres(i,j)))
550	endif !precip
551	!------------------------------------------------------------- endif ! ent/det
552	else !mp nul
553	qDi(i,j,it)=tr(i,j,it) ! orig 0.
554	if(NO_precip(i,j)) then
555	qPr(i,j,it)=0.
556	else
557	qPr(i,j,it)=((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	+imp(i,j)/RGdxpres(i,j)tr(i,j,it))/&
560	(pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))
561	endif
562	endif ! mp nul ou non
563	endif ! condition sur j
564	endif ! kappa
565	ENDDO
566	ENDDO
567
568	!! print test descente insaturee
569	! DO j=klev,1,-1
570	! DO i=1,klon
571	! if(it.eq.3) then
572	! 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,&
573	!! 'zmfdam',zmfdam(i,j,it),'zmfpsi',zmfphi2(i,j,it),&
574	! 'zmfdam+zmfpsi',zmfdam(i,j,it)+zmfphi2(i,j,it),'qpmMint',qpmMint(i,j,it),&
575	! 'Pm',Pm(i,j),'Mint',Mint(i,j),&
576	!! 'zmfa',zmfa(i,j,it),'zmfp',zmfp(i,j,it),&
577	! 'zmfdam',zmfdam(i,j,it),'zmfpsi',zmfphi2(i,j,it),'zmfd1a',zmfd1a(i,j,it)
578	!! 'Pa',Pa(i,j),'eplaMm',eplaMm(i,j)
579	!! 'zmfd1a=da1qa',zmfd1a(i,j,it),'PaqPa',wdtrainA(i,j)*qPa(i,j,it),'da1',d1a(i,j)
580	! endif
581	! ENDDO
582	! ENDDO
583
584
585	! ===================================================
586	! calcul final des tendances
587	! ===================================================
588
589	DO k=klev-1,1,-1
590	DO i=1, klon
591	! transport
592	tdcvMA=zmfd(i,k,it)+zmfu(i,k,it)+zmfa(i,k,it)+zmfp(i,k,it) ! double comptage des downdraft insatures
593	trsptrac=zmfu(i,k,it)+zmfa(i,k,it)+zmfp(i,k,it)
594	! lessivage courants satures
595	scavtrac=-ccntrAA_coef*zmfd1a(i,k,it)&
596	-zmfphi2(i,k,it)*ccntrENV_coef&
597	-zmfdam(i,k,it)*ccntrAA_coef
598	! lessivage courants insatures
599	if(k.le.inb(i).and.k.gt.1) then ! tendances dans le nuage
600	!------------------------------------------------------------- detrainement
601	if(mp(i,k+1).gt.mp(i,k).and.mp(i,k+1).gt.1.e-10) then
602	uscavtrac= (-mp(i,k)+mp(i,k+1))*(qDi(i,k,it)-tr(i,k,it))&
603	+ mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))
604	!
605	! if(it.eq.3) write(*,'(I2,1X,a,5X,e20.12,82X,a,e20.12)')k,' det incloud',&
606	! (-mp(i,k)+mp(i,k+1))(qDi(i,k,it)-tr(i,k,it))pdtimeRG/dxpres(i,k)+&
607	! mp(i,k)(tr(i,k-1,it)-tr(i,k,it))pdtimeRG/dxpres(i,k),&
608	! 'mp',mp(i,k)
609	!------------------------------------------------------------- entrainement
610	elseif(mp(i,k).gt.mp(i,k+1).and.mp(i,k).gt.1.e-10) then
611	uscavtrac= mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))
612	!
613	! 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)
614	!=!------------------------------------------------------------- end ent/det
615	else ! mp(i,k+1)=0. et mp(i,k)=0. pluie directement sur l environnement
616
617	if(NO_precip(i,k)) then
618	uscavtrac=0.
619	! 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)
620	else
621	uscavtrac=-imp(i,k)tr(i,k,it)dxpres(i,k)/RG+evap(i,k)qPr(i,k,it)dxpres(i,k)/RG
622	! 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)
623	endif
624	endif ! mp/det/ent
625	!------------------------------------------------------------- premiere couche
626	elseif(k.eq.1) then ! mp(1)=0.
627	if(mp(i,2).gt.1.e-10) then !detrainement
628	uscavtrac= (-0.+mp(i,2))*(qDi(i,k,it)-tr(i,k,it)) !&
629	! + mp(i,2)*(0.-tr(i,k,it))
630	!
631	! if(it.eq.3) write(*,'(I2,1X,a,e20.12,84X,a,e20.12)')k,' 1 det',&
632	! (-0.+mp(i,2))(qDi(i,k,it)-tr(i,k,it))pdtimeRG/dxpres(i,k)+&
633	! mp(i,2)(0.-tr(i,k,it))pdtimeRG/dxpres(i,k),&
634	! 'mp',mp(i,k)
635	else ! mp(2) = 0 = mp(1) pas de descente insaturee, rien ne se passe s'il ne pleut pas, sinon pluie->env
636	if(NO_precip(i,1)) then
637	uscavtrac=0.
638	else
639	uscavtrac=-imp(i,k)tr(i,k,it)dxpres(i,k)/RG+evap(i,k)qPr(i,k,it)dxpres(i,k)/RG
640	endif
641	! 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)
642	endif
643
644	else ! k > INB au-dessus du nuage
645	uscavtrac=0.
646	endif
647
648	! ===== tendances finales ======
649	trsptd(i,k,it)=trsptrac*pdtimeRG/dxpres(i,k) ! td transport sans eau dans courants satures
650	dtrSscav(i,k,it)=scavtrac*pdtimeRG/dxpres(i,k) ! td du lessivage dans courants satures
651	dtrUscav(i,k,it)=uscavtrac*pdtimeRG/dxpres(i,k) ! td courant insat
652	dtrsat(i,k,it)=(trsptrac+scavtrac)*pdtimeRG/dxpres(i,k) ! td courant sat
653	dtrcv(i,k,it)=(trsptrac+scavtrac+uscavtrac)*pdtimeRG/dxpres(i,k)!dtrsat(i,k,it)+dtrUscav(i,k,it) td conv
654	!!!!!!
655	dtrcvMA(i,k,it)=tdcvMA*pdtimeRG/dxpres(i,k) ! MA tendance convection
656	ENDDO
657	ENDDO
658
659	! test de conservation du traceur
660	!print*,"_____________________________________________________________"
661	!print*," "
662	! conserv=0.
663	! conservMA=0.
664	! DO k= klev-1,1,-1
665	! DO i=1, klon
666	! conserv=conserv+dtrcv(i,k,it)* &
667	! (paprs(i,k)-paprs(i,k+1))/RG
668	! conservMA=conservMA+dtrcvMA(i,k,it)* &
669	! (paprs(i,k)-paprs(i,k+1))/RG
670	!
671	! if(it.eq.3) write(*,'(I2,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12)') k,&
672	! 'MA td ',dtrcvMA(i,k,it)*dxpres(i,k)/RG,&
673	! ' td',dtrcv(i,k,it)*dxpres(i,k)/RG,' conservMA ',conservMA,'conserv ',conserv
674	!!
675	! ENDDO
676	! ENDDO
677	! if(it.eq.3) print *,'it',it,'conserv ',conserv,'conservMA ',conservMA
678
679	END SUBROUTINE cvltr

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

Download in other formats:

Original Format