Context Navigation

cvltr.f90 @ 5441

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

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

Download in other formats:

Original Format