Context Navigation

cvltr_spl.f90 @ 5301

Last change on this file since 5301 was 5292, checked in by abarral, 4 days ago
Move academic.h chem.h chem_spla.h to module
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 chem_mod_h
14	USE yoecumf_mod_h
15	USE conema3_mod_h
16	USE IOIPSL
17	USE dimphy
18	USE infotrac_phy, ONLY : nbtr
19	USE yomcst_mod_h
20	IMPLICIT NONE
21	!=====================================================================
22	! Objet : convection des traceurs / KE
23	! Auteurs: M-A Filiberti and J-Y Grandpeix
24	! modifiee par R Pilon : lessivage des traceurs / KE
25	!=====================================================================
26
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	!
215	! initialisation pour flux de traceurs, td et autre
216	trsptrac = 0.
217	scavtrac = 0.
218	uscavtrac = 0.
219	qfeed(:,it) = 0. !RL
220	DO j=1,klev
221	DO i=1,klon
222	zmfd(i,j,it)=0.
223	zmfa(i,j,it)=0.
224	zmfu(i,j,it)=0.
225	zmfp(i,j,it)=0.
226	zmfphi2(i,j,it)=0.
227	zmfd1a(i,j,it)=0.
228	zmfdam(i,j,it)=0.
229	qDi(i,j,it)=0.
230	qPr(i,j,it)=0.
231	qPa(i,j,it)=0.
232	qMel(i,j,it)=0.
233	qMeltmp(i,j,it)=0.
234	qTrdi(i,j,it)=0.
235	kappa(i,j)=0.
236	trsptd(i,j,it)=0.
237	dtrsat(i,j,it)=0.
238	dtrSscav(i,j,it)=0.
239	dtrUscav(i,j,it)=0.
240	dtrcv(i,j,it)=0.
241	dtrcvMA(i,j,it)=0.
242	evap(i,j)=0.
243	dxpres(i,j)=0.
244	qpmMint(i,j,it)=0.
245	Mint(i,j)=0.
246	END DO
247	END DO
248
249	! suppression des valeurs très faibles (~1e-320)
250	! multiplication de levaporation pour lavoir par unite de temps
251	! et par unite de surface de la maille
252	! -> cv30_unsat : evap : masse evaporee/s/(m2 de la descente)
253	DO j=1,klev
254	DO i=1,klon
255	if(ev(i,j).lt.1.e-16) then
256	evap(i,j)=0.
257	else
258	evap(i,j)=ev(i,j)*sigd(i)
259	endif
260	END DO
261	END DO
262
263	DO j=1,klev
264	DO i=1,klon
265	if(j.lt.klev) then
266	if(epIN(i,j).lt.1.e-32) then
267	ep(i,j)=0.
268	else
269	ep(i,j)=epIN(i,j)
270	endif
271	else
272	ep(i,j)=epmax
273	endif
274	if(mpIN(i,j).lt.1.e-32) then
275	mp(i,j)=0.
276	else
277	mp(i,j)=mpIN(i,j)
278	endif
279	if(pmflxsIN(i,j).lt.1.e-32) then
280	pmflxs(i,j)=0.
281	else
282	pmflxs(i,j)=pmflxsIN(i,j)
283	endif
284	if(pmflxrIN(i,j).lt.1.e-32) then
285	pmflxr(i,j)=0.
286	else
287	pmflxr(i,j)=pmflxrIN(i,j)
288	endif
289	if(wdtrainA(i,j).lt.1.e-32) then
290	Pa(i,j)=0.
291	else
292	Pa(i,j)=wdtrainA(i,j)
293	endif
294	if(wdtrainM(i,j).lt.1.e-32) then
295	Pm(i,j)=0.
296	else
297	Pm(i,j)=wdtrainM(i,j)
298	endif
299	END DO
300	END DO
301
302	!==========================================
303	DO j = klev-1,1,-1
304	DO i = 1,klon
305	NO_precip(i,j) = (pmflxr(i,j+1)+pmflxs(i,j+1)).lt.1.e-10&
306	.and.Pa(i,j).lt.1.e-10.and.Pm(i,j).lt.1.e-10
307	END DO
308	END DO
309	!==============================================================================
310	! JE SPLA: Calc of ccntrAA_coef,ccntrENV_coef, coefcoli for SPLA aerosols and
311	! precursors. From SPLA code inscav_spl.F
312	!print *,'Using SPLA values for cvltr_spl ccntr and coefcoli params'
313	!
314	!
315	!IF (it.EQ.2) THEN !--fine aerosol
316	! DO j=1,klev
317	! DO i=1,klon
318	! ccntrAA_coef3d(i,j)=0.7
319	! ccntrENV_coef3d(i,j)=0.7
320	! coefcoli3d(i,j)=0.001
321	! ENDDO
322	! ENDDO
323	!ELSEIF (it.EQ.3) THEN !-- Coarse Sea salt aerosol
324	! DO j=1,klev
325	! DO i=1,klon
326	! ccntrAA_coef3d(i,j)=1.0
327	! ccntrENV_coef3d(i,j)=1.0
328	! coefcoli3d(i,j)=0.001
329	! ENDDO
330	! ENDDO
331	!
332	!ELSEIF (it.EQ.4) THEN !--Coarse Dust aerosol
333	! DO j=1,klev
334	! DO i=1,klon
335	! ccntrAA_coef3d(i,j)=0.7
336	! ccntrENV_coef3d(i,j)=0.7
337	! coefcoli3d(i,j)=0.001
338	!
339	! ENDDO
340	! ENDDO
341	! Gas precursor: Henry's law
342
343	IF (it .EQ. id_prec) THEN
344	DO k=1, klev
345	DO i=1, klon
346	henry_t=henry(it)exp(-kk(it)(1./298.-1./te(i,k))) !--mol/l/atm
347	K1=1.2e-2exp(-2010(1/298.-1/te(i,k)))
348	K2=6.6e-8exp(-1510(1/298.-1/te(i,k)))
349	henry_t=henry_t(1. + K1/10.(-ph) + K1K2/(10.(-ph))2)
350	f_a=henry_t/101.325Rte(i,k)qliqzrho(i,k)/rho_water
351	! scav(i,k)=f_a/(1.+f_a)
352	ccntrAA_coef3d(i,k)= f_a/(1.+f_a)
353	ccntrENV_coef3d(i,k)= f_a/(1.+f_a)
354	coefcoli3d(i,k)=0.0
355	ENDDO
356	ENDDO
357	! CALL minmaxqfi2(clw,1.e33,-1.e33,'minmax clw')
358	ELSE
359	DO j=1,klev
360	DO i=1,klon
361	ccntrAA_coef3d(i,j)=ccntrAA_spla(it)
362	ccntrENV_coef3d(i,j)=ccntrENV_spla(it)
363	coefcoli3d(i,j)=coefcoli_spla(it)
364	ENDDO
365	ENDDO
366
367
368	ENDIF
369
370	! JE end SPLA modifs in ccn parameters
371	!==============================================================================
372
373
374
375
376	!JE init 20140103
377	! impaction sur la surface de la colonne de la descente insaturee
378	! On prend la moyenne des precip entre le niveau i+1 et i
379	! I=3/4* (P(1+1)+P(i))/2 / (sigdrrho_l)
380	! 1000kg/m3= densité de l'eau
381	! 0.75e-3 = 3/4 /1000
382	! Par la suite, I est tout le temps multiplié par sig_d pour avoir l'impaction sur la surface de la maille
383	! on le néglige ici pour simplifier le code
384	do j=1,klev-1
385	do i=1,klon
386	!JE imp(i,j) = coefcoli0.75e-3/rdrop &
387	!JE 0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j))
388	imp(i,j) = coefcoli3d(i,j)0.75e-3/rdrop &
389	0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j))
390	! rho(i,j)=pplay(i,j)/(rd*te(i,j))
391	enddo
392	enddo
393	!JE end 20140103
394
395
396	! =========================================
397	! calcul des tendances liees au downdraft
398	! =========================================
399	!cdir collapse
400	DO k=1,klev
401	DO j=1,klev
402	DO i=1,klon
403	zmd(i,j,k)=0.
404	za (i,j,k)=0.
405	END DO
406	END DO
407	END DO
408	! calcul de la matrice d echange
409	! matrice de distribution de la masse entrainee en k
410	! commmentaire RomP : mp > 0
411	DO k=1,klev-1
412	DO i=1,klon
413	zmd(i,k,k)=max(0.,mp(i,k)-mp(i,k+1)) ! ~ mk(k)
414	END DO
415	END DO
416	DO k=2,klev
417	DO j=k-1,1,-1
418	DO i=1,klon
419	if(mp(i,j+1).gt.1.e-10) then
420	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)
421	ENDif
422	END DO
423	END DO
424	END DO
425	DO k=1,klev
426	DO j=1,klev-1
427	DO i=1,klon
428	za(i,j,k)=max(0.,zmd(i,j+1,k)-zmd(i,j,k))
429	END DO
430	END DO
431	END DO
432	!!!!! quantite de traceur dans la descente d'air insaturee : 4 juin 2012
433	DO k=1,klev
434	DO j=1,klev-1
435	DO i=1,klon
436	if(mp(i,j+1).gt.1.e-10) then
437	qTrdi(i,j+1,it)=qTrdi(i,j+1,it)+(zmd(i,j+1,k)/mp(i,j+1))*tr(i,k,it)
438	else
439	qTrdi(i,j,it)=0.!tr(i,j,it)
440	endif
441	ENDDO
442	ENDDO
443	ENDDO
444	!!!!!
445	!
446	! rajout du terme lie a l ascendance induite
447	!
448	DO j=2,klev
449	DO i=1,klon
450	za(i,j,j-1)=za(i,j,j-1)+mp(i,j)
451	END DO
452	END DO
453	!
454	! tendance courants insatures ! sans lessivage ancien schema
455	!
456	DO k=1,klev
457	DO j=1,klev
458	DO i=1,klon
459	zmfd(i,j,it)=zmfd(i,j,it)+za(i,j,k)*(tr(i,k,it)-tr(i,j,it))
460	END DO
461	END DO
462	END DO
463	!
464	! =========================================
465	! calcul des tendances liees aux courants satures j <-> z ; k <-> z'
466	! =========================================
467	!RL
468	! Feeding concentrations
469	DO j=1,klev
470	DO i=1,klon
471	qfeed(i,it)=qfeed(i,it)+wght_cvfd(i,j)*tr(i,j,it)
472	END DO
473	END DO
474	!RL
475	!
476	DO j=1,klev
477	DO i=1,klon
478	!RL
479	!! zmfa(i,j,it)=da(i,j)*(tr(i,1,it)-tr(i,j,it)) ! da
480	zmfa(i,j,it)=da(i,j)*(qfeed(i,it)-tr(i,j,it)) ! da
481	!RL
482	END DO
483	END DO
484	!
485	DO k=1,klev
486	DO j=1,klev
487	DO i=1,klon
488	zmfp(i,j,it)=zmfp(i,j,it)+phi(i,j,k)*(tr(i,k,it)-tr(i,j,it)) ! phi
489	END DO
490	END DO
491	END DO
492	! RomP ajout des matrices liees au lessivage
493	DO j=1,klev
494	DO i=1,klon
495	zmfd1a(i,j,it)=d1a(i,j)*tr(i,1,it) ! da1
496	zmfdam(i,j,it)=dam(i,j)*tr(i,1,it) ! dam
497	END DO
498	END DO
499	DO k=1,klev
500	DO j=1,klev
501	DO i=1,klon
502	zmfphi2(i,j,it)=zmfphi2(i,j,it)+phi2(i,j,k)*tr(i,k,it) ! psi
503	END DO
504	END DO
505	END DO
506	DO j=1,klev-1
507	DO i=1,klon
508	zmfu(i,j,it)=max(0.,upd(i,j+1)+dnd(i,j+1))*(tr(i,j+1,it)-tr(i,j,it))
509	END DO
510	END DO
511	DO j=2,klev
512	DO i=1,klon
513	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))
514	END DO
515	END DO
516	! ===================================================
517	! calcul des tendances liees aux courants insatures
518	! ===================================================
519	! pression
520	DO k=1, klev
521	DO i=1, klon
522	dxpres(i,k)=paprs(i,k)-paprs(i,k+1)
523	ENDDO
524	ENDDO
525	pdtimeRG=pdtime*RG
526
527	! q_pa et q_pm traceurs issues des courants satures se retrouvant dans les precipitations
528	DO j=1,klev
529	DO i=1,klon
530	if(j.ge.icb(i).and.j.le.inb(i)) then
531	if(clw(i,j).gt.1.e-16) then
532	! qPa(i,j,it)=ccntrAA_coef*tr(i,1,it)/clw(i,j)
533	qPa(i,j,it)=ccntrAA_coef3d(i,j)*tr(i,1,it)/clw(i,j)
534	else
535	qPa(i,j,it)=0.
536	endif
537	endif
538	END DO
539	END DO
540
541	! calcul de q_pm en 2 parties :
542	! 1) calcul de sa valeur pour un niveau z' donne
543	! 2) integration sur la verticale sur z'
544	DO j=1,klev
545	DO k=1,j-1
546	DO i=1,klon
547	if(k.ge.icb(i).and.k.le.inb(i).and.&
548	j.le.inb(i)) then
549	if(elij(i,k,j).gt.1.e-16) then
550	!JE qMeltmp(i,j,it)=((1-ep(i,k))ccntrAA_coeftr(i,1,it)&
551	!JE *(1.-sij(i,k,j)) +ccntrENV_coef&
552	!JE tr(i,k,it)sij(i,k,j)) / elij(i,k,j)
553	qMeltmp(i,j,it)=((1-ep(i,k))ccntrAA_coef3d(i,k)tr(i,1,it)&
554	*(1.-sij(i,k,j)) +ccntrENV_coef3d(i,k)&
555	tr(i,k,it)sij(i,k,j)) / elij(i,k,j)
556	else
557	qMeltmp(i,j,it)=0.
558	endif
559	qpmMint(i,j,it)=qpmMint(i,j,it) + qMeltmp(i,j,it)*epmlmMm(i,j,k)
560	Mint(i,j)=Mint(i,j) + epmlmMm(i,j,k)
561	endif ! end if dans nuage
562	END DO
563	END DO
564	END DO
565
566	DO j=1,klev
567	DO i=1,klon
568	if(Mint(i,j).gt.1.e-16) then
569	qMel(i,j,it)=qpmMint(i,j,it)/Mint(i,j)
570	else
571	qMel(i,j,it)=0.
572	endif
573	END DO
574	END DO
575
576	! calcul de q_d et q_p traceurs de la descente precipitante
577	DO j=klev-1,1,-1
578	DO i=1,klon
579	if(mp(i,j+1).gt.mp(i,j).and.mp(i,j+1).gt.1.e-10) then ! detrainement
580	kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
581	(-mp(i,j+1)-imp(i,j)/RG*dxpres(i,j))&
582	+ (imp(i,j)/RGdxpres(i,j))(evap(i,j)/RG*dxpres(i,j)))
583
584	elseif(mp(i,j).gt.mp(i,j+1).and.mp(i,j).gt.1.e-10) then! entrainement
585	if(j.eq.1) then
586	kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
587	(-mp(i,2)-imp(i,j)/RG*dxpres(i,j))&
588	+ (imp(i,j)/RGdxpres(i,j))(evap(i,j)/RG*dxpres(i,j)))
589	else
590	kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
591	(-mp(i,j)-imp(i,j)/RG*dxpres(i,j))&
592	+ (imp(i,j)/RGdxpres(i,j))(evap(i,j)/RG*dxpres(i,j)))
593	endif
594	else
595	kappa(i,j)=1.
596	endif
597	ENDDO
598	ENDDO
599
600	DO j=klev-1,1,-1
601	DO i=1,klon
602	if (abs(kappa(i,j)).lt.1.e-25) then !si denominateur nul (il peut y avoir des mp!=0)
603	kappa(i,j)=1.
604	if(j.eq.1) then
605	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
606	elseif(mp(i,j+1).gt.mp(i,j).and.mp(i,j+1).gt.1.e-10) then
607	qDi(i,j,it)=qDi(i,j+1,it)
608	elseif(mp(i,j).gt.mp(i,j+1).and.mp(i,j).gt.1.e-10) then! entrainement
609	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))
610	else ! si mp (i)=0 et mp(j+1)=0
611	qDi(i,j,it)=tr(i,j,it) ! orig 0.
612	endif
613
614	if(NO_precip(i,j)) then
615	qPr(i,j,it)=0.
616	else
617	qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
618	Pa(i,j)qPa(i,j,it)+Pm(i,j)qMel(i,j,it)&
619	+imp(i,j)/RGdxpres(i,j)qDi(i,j,it))/&
620	(pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))
621	endif
622	else ! denominateur non nul
623	kappa(i,j)=1./kappa(i,j)
624	! calcul de qd et qp
625	!!jyg (20130119) correction pour le sommet du nuage
626	!! if(j.ge.inb(i)) then !au-dessus du nuage, sommet inclu
627	if(j.gt.inb(i)) then !au-dessus du nuage
628	qDi(i,j,it)=tr(i,j,it) ! pas de descente => environnement = descente insaturee
629	qPr(i,j,it)=0.
630
631	! vvv premiere couche du modele ou mp(1)=0 ! det tout le temps vvv
632	elseif(j.eq.1) then
633	if(mp(i,2).gt.1.e-10) then !mp(2) non nul -> detrainement (car mp(1) = 0) !ent pas possible
634	if(NO_precip(i,j)) then !pas de precip en (i)
635	qDi(i,j,it)=qDi(i,j+1,it)
636	qPr(i,j,it)=0.
637	else
638	qDi(i,j,it)=kappa(i,j)*(&
639	(-evap(i,j)/RGdxpres(i,j))((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
640	Pa(i,j)qPa(i,j,it)+Pm(i,j)qMel(i,j,it)) +&
641	(pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
642	(-mp(i,j+1)*qDi(i,j+1,it)))
643
644	qPr(i,j,it)=kappa(i,j)*(&
645	(-mp(i,j+1)-imp(i,j)/RGdxpres(i,j))&
646	((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
647	Pa(i,j)qPa(i,j,it)+Pm(i,j)qMel(i,j,it))&
648	+(-mp(i,j+1)qDi(i,j+1,it)) (imp(i,j)/RG*dxpres(i,j)))
649	endif
650
651	else !mp(2) nul -> plus de descente insaturee -> pluie agit sur environnement
652	qDi(i,j,it)=tr(i,j,it) ! orig 0.
653	if(NO_precip(i,j)) then
654	qPr(i,j,it)=0.
655	else
656	qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
657	Pa(i,j)qPa(i,j,it)+Pm(i,j)qMel(i,j,it)&
658	+imp(i,j)/RGdxpres(i,j)tr(i,j,it))/&
659	(pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))
660	endif
661
662	endif !mp(2) nul ou non
663
664	! vvv (j!=1.and.j.lt.inb(i)) en-dessous du sommet nuage vvv
665	else
666	!------------------------------------------------------------- detrainement
667	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
668	if(NO_precip(i,j)) then
669	qDi(i,j,it)=qDi(i,j+1,it)
670	qPr(i,j,it)=0.
671	else
672	qDi(i,j,it)=kappa(i,j)*(&
673	(-evap(i,j)/RGdxpres(i,j))((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
674	Pa(i,j)qPa(i,j,it)+Pm(i,j)qMel(i,j,it)) +&
675	(pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
676	(-mp(i,j+1)*qDi(i,j+1,it)))
677	!
678	qPr(i,j,it)=kappa(i,j)*(&
679	(-mp(i,j+1)-imp(i,j)/RGdxpres(i,j))&
680	((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
681	Pa(i,j)qPa(i,j,it)+Pm(i,j)qMel(i,j,it))&
682	+(-mp(i,j+1)qDi(i,j+1,it)) (imp(i,j)/RG*dxpres(i,j)))
683	endif !precip
684	!------------------------------------------------------------- entrainement
685	elseif(mp(i,j).gt.mp(i,j+1).and.mp(i,j).gt.1.e-10) then
686	if(NO_precip(i,j)) then
687	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))
688	qPr(i,j,it)=0.
689	else
690	qDi(i,j,it)=kappa(i,j)*(&
691	(-evap(i,j)/RGdxpres(i,j))((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
692	Pa(i,j)qPa(i,j,it)+Pm(i,j)qMel(i,j,it)) +&
693	(pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*&
694	(-mp(i,j+1)(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)tr(i,j,it)))
695	!
696	qPr(i,j,it)=kappa(i,j)*(&
697	(-mp(i,j)-imp(i,j)/RGdxpres(i,j))&
698	((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
699	Pa(i,j)qPa(i,j,it)+Pm(i,j)qMel(i,j,it))&
700	+(-mp(i,j+1)(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)tr(i,j,it))*&
701	(imp(i,j)/RG*dxpres(i,j)))
702	endif !precip
703	!------------------------------------------------------------- endif ! ent/det
704	else !mp nul
705	qDi(i,j,it)=tr(i,j,it) ! orig 0.
706	if(NO_precip(i,j)) then
707	qPr(i,j,it)=0.
708	else
709	qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+&
710	Pa(i,j)qPa(i,j,it)+Pm(i,j)qMel(i,j,it)&
711	+imp(i,j)/RGdxpres(i,j)tr(i,j,it))/&
712	(pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))
713	endif
714	endif ! mp nul ou non
715	endif ! condition sur j
716	endif ! kappa
717	ENDDO
718	ENDDO
719
720	!! print test descente insaturee
721	! DO j=klev,1,-1
722	! DO i=1,klon
723	! if(it.eq.3) then
724	! 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,&
725	!! 'zmfdam',zmfdam(i,j,it),'zmfpsi',zmfphi2(i,j,it),&
726	! 'zmfdam+zmfpsi',zmfdam(i,j,it)+zmfphi2(i,j,it),'qpmMint',qpmMint(i,j,it),&
727	! 'Pm',Pm(i,j),'Mint',Mint(i,j),&
728	!! 'zmfa',zmfa(i,j,it),'zmfp',zmfp(i,j,it),&
729	! 'zmfdam',zmfdam(i,j,it),'zmfpsi',zmfphi2(i,j,it),'zmfd1a',zmfd1a(i,j,it)
730	!! 'Pa',Pa(i,j),'eplaMm',eplaMm(i,j)
731	!! 'zmfd1a=da1qa',zmfd1a(i,j,it),'PaqPa',wdtrainA(i,j)*qPa(i,j,it),'da1',d1a(i,j)
732	! endif
733	! ENDDO
734	! ENDDO
735
736
737	! ===================================================
738	! calcul final des tendances
739	! ===================================================
740
741	DO k=klev-1,1,-1
742	DO i=1, klon
743	! transport
744	tdcvMA=zmfd(i,k,it)+zmfu(i,k,it)+zmfa(i,k,it)+zmfp(i,k,it) ! double comptage des downdraft insatures
745	trsptrac=zmfu(i,k,it)+zmfa(i,k,it)+zmfp(i,k,it)
746	! lessivage courants satures
747	!JE scavtrac=-ccntrAA_coef*zmfd1a(i,k,it)&
748	!JE -zmfphi2(i,k,it)*ccntrENV_coef&
749	!JE -zmfdam(i,k,it)*ccntrAA_coef
750	scavtrac=-ccntrAA_coef3d(i,k)*zmfd1a(i,k,it)&
751	-zmfphi2(i,k,it)*ccntrENV_coef3d(i,k)&
752	-zmfdam(i,k,it)*ccntrAA_coef3d(i,k)
753	! lessivage courants insatures
754	if(k.le.inb(i).and.k.gt.1) then ! tendances dans le nuage
755	!------------------------------------------------------------- detrainement
756	if(mp(i,k+1).gt.mp(i,k).and.mp(i,k+1).gt.1.e-10) then
757	uscavtrac= (-mp(i,k)+mp(i,k+1))*(qDi(i,k,it)-tr(i,k,it))&
758	+ mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))
759	!
760	! if(it.eq.3) write(*,'(I2,1X,a,5X,e20.12,82X,a,e20.12)')k,' det incloud',&
761	! (-mp(i,k)+mp(i,k+1))(qDi(i,k,it)-tr(i,k,it))pdtimeRG/dxpres(i,k)+&
762	! mp(i,k)(tr(i,k-1,it)-tr(i,k,it))pdtimeRG/dxpres(i,k),&
763	! 'mp',mp(i,k)
764	!------------------------------------------------------------- entrainement
765	elseif(mp(i,k).gt.mp(i,k+1).and.mp(i,k).gt.1.e-10) then
766	uscavtrac= mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))
767	!
768	! 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)
769	!=!------------------------------------------------------------- end ent/det
770	else ! mp(i,k+1)=0. et mp(i,k)=0. pluie directement sur l environnement
771
772	if(NO_precip(i,k)) then
773	uscavtrac=0.
774	! 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)
775	else
776	uscavtrac=-imp(i,k)tr(i,k,it)dxpres(i,k)/RG+evap(i,k)qPr(i,k,it)dxpres(i,k)/RG
777	! 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)
778	!!JE adds
779	! if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'imp',imp(i,k)
780	! 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)
781	! if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'evap',evap(i,k)
782	! 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)
783	! if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'dxpres',dxpres(i,k)
784	!!Je end
785
786	endif
787	endif ! mp/det/ent
788	!------------------------------------------------------------- premiere couche
789	elseif(k.eq.1) then ! mp(1)=0.
790	if(mp(i,2).gt.1.e-10) then !detrainement
791	uscavtrac= (-0.+mp(i,2))*(qDi(i,k,it)-tr(i,k,it)) !&
792	! + mp(i,2)*(0.-tr(i,k,it))
793	!
794	! if(it.eq.3) write(*,'(I2,1X,a,e20.12,84X,a,e20.12)')k,' 1 det',&
795	! (-0.+mp(i,2))(qDi(i,k,it)-tr(i,k,it))pdtimeRG/dxpres(i,k)+&
796	! mp(i,2)(0.-tr(i,k,it))pdtimeRG/dxpres(i,k),&
797	! 'mp',mp(i,k)
798	else ! mp(2) = 0 = mp(1) pas de descente insaturee, rien ne se passe s'il ne pleut pas, sinon pluie->env
799	if(NO_precip(i,1)) then
800	uscavtrac=0.
801	else
802	uscavtrac=-imp(i,k)tr(i,k,it)dxpres(i,k)/RG+evap(i,k)qPr(i,k,it)dxpres(i,k)/RG
803	endif
804	! 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)
805	endif
806
807	else ! k > INB au-dessus du nuage
808	uscavtrac=0.
809	endif
810
811	! ===== tendances finales ======
812	trsptd(i,k,it)=trsptrac*pdtimeRG/dxpres(i,k) ! td transport sans eau dans courants satures
813	dtrSscav(i,k,it)=scavtrac*pdtimeRG/dxpres(i,k) ! td du lessivage dans courants satures
814	dtrUscav(i,k,it)=uscavtrac*pdtimeRG/dxpres(i,k) ! td courant insat
815	dtrsat(i,k,it)=(trsptrac+scavtrac)*pdtimeRG/dxpres(i,k) ! td courant sat
816	dtrcv(i,k,it)=(trsptrac+scavtrac+uscavtrac)*pdtimeRG/dxpres(i,k)!dtrsat(i,k,it)+dtrUscav(i,k,it) td conv
817	!!!!!!
818	dtrcvMA(i,k,it)=tdcvMA*pdtimeRG/dxpres(i,k) ! MA tendance convection
819	ENDDO
820	ENDDO
821
822	! test de conservation du traceur
823	!print*,"_____________________________________________________________"
824	!print*," "
825	! conserv=0.
826	! conservMA=0.
827	! DO k= klev-1,1,-1
828	! DO i=1, klon
829	! conserv=conserv+dtrcv(i,k,it)* &
830	! (paprs(i,k)-paprs(i,k+1))/RG
831	! conservMA=conservMA+dtrcvMA(i,k,it)* &
832	! (paprs(i,k)-paprs(i,k+1))/RG
833	!
834	! if(it.eq.3) write(*,'(I2,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12)') k,&
835	! 'MA td ',dtrcvMA(i,k,it)*dxpres(i,k)/RG,&
836	! ' td',dtrcv(i,k,it)*dxpres(i,k)/RG,' conservMA ',conservMA,'conserv ',conserv
837	!!
838	! ENDDO
839	! ENDDO
840	! if(it.eq.3) print *,'it',it,'conserv ',conserv,'conservMA ',conservMA
841
842	END SUBROUTINE cvltr_spl

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

Download in other formats:

Original Format