Context Navigation

cvltr.F90 @ 5371

Last change on this file since 5371 was 2320, checked in by Ehouarn Millour, 9 years ago
Physics/dynamics separation: make an infotrac_phy module, which should be used from within the physics, and is initialized from infotrac (dynamics) via iniphysiq. EM
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.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