Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

cvltr_spl.F90

Last change on this file was 5160, checked in by abarral, 3 months ago
Put .h into modules
File size: 30.2 KB

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

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

Download in other formats:

Original Format