Context Navigation

cvltr.F90 @ 2942

Last change on this file since 2942 was 1759, checked in by Ehouarn Millour, 11 years ago
IOIPSL routine getin is not threadsafe, so when running in OpenMP, it should be called by only one thread (and the result copied to other threads in the case of threadprivate variables). EM
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