Context Navigation

cvltr.F90 @ 3814

Last change on this file since 3814 was 3814, checked in by ymipsl, 10 years ago

remove all dynamic dependency in LMDZ physics except for the include "dimensions.h"

File size: 25.6 KB

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

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

Download in other formats:

Original Format