Context Navigation

cvltr_spl.F90 @ 3627

Last change on this file since 3627 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
File size: 30.3 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,tname
16	IMPLICIT NONE
17	!=====================================================================
18	! Objet : convection des traceurs / KE
19	! Auteurs: M-A Filiberti and J-Y Grandpeix
20	! modifiee par R Pilon : lessivage des traceurs / KE
21	!=====================================================================
22
23	include "YOMCST.h"
24	include "YOECUMF.h"
25	include "conema3.h"
26	include "chem.h"
27
28	! Entree
29	REAL,INTENT(IN) :: pdtime
30	REAL,DIMENSION(klon,klev),INTENT(IN) :: da
31	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi
32	! RomP
33	REAL,DIMENSION(klon,klev),INTENT(IN) :: d1a,dam ! matrices pour simplifier
34	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi2 ! l'ecriture des tendances
35	!
36	REAL,DIMENSION(klon,klev),INTENT(IN) :: mpIN
37	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression aux 1/2 couches (bas en haut)
38	! REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression aux 1/2 couches (bas en haut)
39	REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr ! q de traceur (bas en haut)
40	INTEGER,INTENT(IN) :: it
41	REAL,DIMENSION(klon,klev),INTENT(IN) :: upd ! saturated updraft mass flux
42	REAL,DIMENSION(klon,klev),INTENT(IN) :: dnd ! saturated downdraft mass flux
43	!
44	REAL,DIMENSION(klon,klev),INTENT(IN) :: wdtrainA ! masses precipitantes de l'asc adiab
45	REAL,DIMENSION(klon,klev),INTENT(IN) :: wdtrainM ! masses precipitantes des melanges
46	!JE REAL,DIMENSION(klon,klev),INTENT(IN) :: pmflxrIN ! vprecip: eau
47	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: pmflxrIN ! vprecip: eau
48	!JE REAL,DIMENSION(klon,klev),INTENT(IN) :: pmflxsIN ! vprecip: neige
49	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: pmflxsIN ! vprecip: neige
50	REAL,DIMENSION(klon,klev),INTENT(IN) :: ev ! evaporation cv30_routine
51	REAL,DIMENSION(klon,klev),INTENT(IN) :: epIN
52	REAL,DIMENSION(klon,klev),INTENT(IN) :: te
53	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: sij ! fraction dair de lenv
54	REAL,DIMENSION(klon,klev),INTENT(IN) :: wght_cvfd ! weights of the layers feeding convection
55	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: elij ! contenu en eau condensée spécifique/conc deau condensée massique
56	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: epmlmMm ! eau condensee precipitee dans mel masse dair sat
57	REAL,DIMENSION(klon,klev),INTENT(IN) :: eplaMm ! eau condensee precipitee dans aa masse dair sat
58
59	REAL,DIMENSION(klon,klev),INTENT(IN) :: clw ! contenu en eau condensée dans lasc adiab
60	REAL,DIMENSION(klon),INTENT(IN) :: sigd
61	INTEGER,DIMENSION(klon),INTENT(IN) :: icb,inb
62	! Sortie
63	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrcv ! tendance totale (bas en haut)
64	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrcvMA ! M-A Filiberti
65	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: trsptd ! tendance du transport
66	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrSscav ! tendance du lessivage courant sat
67	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrsat ! tendance trsp+sat scav
68	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrUscav ! tendance du lessivage courant unsat
69	!
70	! Variables locales
71	INTEGER :: i,j,k
72	REAL,DIMENSION(klon,klev) :: dxpres ! difference de pression entre niveau (j+1) et (j)
73	REAL :: pdtimeRG ! pas de temps * gravite
74	REAL,DIMENSION(klon,nbtr) :: qfeed ! tracer concentration feeding convection
75	! variables pour les courants satures
76	REAL,DIMENSION(klon,klev,klev) :: zmd
77	REAL,DIMENSION(klon,klev,klev) :: za
78	REAL,DIMENSION(klon,klev,nbtr) :: zmfd,zmfa
79	REAL,DIMENSION(klon,klev,nbtr) :: zmfp,zmfu
80
81	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfd1a
82	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfdam
83	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfphi2
84
85	! RomP ! les variables sont nettoyees des valeurs aberrantes
86	REAL,DIMENSION(klon,klev) :: Pa, Pm ! pluie AA et mélanges, var temporaire
87	REAL,DIMENSION(klon,klev) :: pmflxs,pmflxr ! pmflxrIN,pmflxsIN sans valeur aberante
88	REAL,DIMENSION(klon,klev) :: mp ! flux de masse
89	REAL,DIMENSION(klon,klev) :: ep ! fraction d'eau convertie en precipitation
90	REAL,DIMENSION(klon,klev) :: evap ! evaporation : variable temporaire
91	REAL,DIMENSION(klon,klev) :: rho !environmental density
92
93	REAL,DIMENSION(klon,klev) :: kappa ! denominateur du au calcul de la matrice
94	! pour obtenir qd et qp
95	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qTrdi ! traceurs descente air insature transport MA
96	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qDi ! traceurs descente insaturees
97	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qPr ! traceurs colonne precipitante
98	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qPa ! traceurs dans les precip issues lasc. adiab.
99	REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qMel ! traceurs dans les precip issues des melanges
100	REAL,DIMENSION(klon,klev,nbtr) :: qMeltmp ! variable temporaire
101	REAL,DIMENSION(klon,klev,nbtr) :: qpmMint
102	REAL,DIMENSION(klon,klev),INTENT(OUT) :: Mint
103	! tendances
104	REAL :: tdcvMA ! terme de transport de traceur (schema Marie Angele)
105	REAL :: trsptrac ! terme de transport de traceur par l'air
106	REAL :: scavtrac ! terme de lessivage courant sature
107	REAL :: uscavtrac ! terme de lessivage courant insature
108	! impaction
109	!!! Correction apres discussion Romain P. / Olivier B.
110	!!! REAL,PARAMETER :: rdrop=2.5e-3 ! rayon des gouttes d'eau
111	REAL,PARAMETER :: rdrop=1.e-3 ! rayon des gouttes d'eau
112	!!!
113	REAL,DIMENSION(klon,klev) :: imp ! coefficient d'impaction
114	! parametres lessivage
115	REAL :: ccntrAA_coef ! \alpha_a : fract aerosols de l'AA convertis en CCN
116	REAL :: ccntrENV_coef ! \beta_m : fract aerosols de l'env convertis en CCN
117	REAL :: coefcoli ! coefficient de collision des gouttes sur les aerosols
118	!
119	LOGICAL,DIMENSION(klon,klev) :: NO_precip
120	! LOGICAL :: scavON
121	! var tmp tests
122	REAL :: conserv
123	real :: conservMA
124	! JE SPLA adaptation
125
126	INTEGER :: id_prec,id_fine,id_coss, id_codu, id_scdu
127	REAL,DIMENSION(nbtr) :: ccntrAA_spla,ccntrENV_spla,coefcoli_spla
128	REAL,DIMENSION(klon,klev) :: ccntrAA_coef3d
129	REAL,DIMENSION(klon,klev) :: ccntrENV_coef3d
130	REAL,DIMENSION(klon,klev) :: coefcoli3d
131
132	REAL,DIMENSION(nbtr) :: henry !--cste de Henry mol/l/atm
133	REAL,DIMENSION(nbtr) :: kk !--coefficient de var avec T (K)
134	REAL :: henry_t !--constante de Henry a T t (mol/l/atm)
135	REAL :: f_a !--rapport de la phase aqueuse a la phase gazeuse
136
137	REAL, PARAMETER :: ph=5.
138	REAL :: K1, K2
139	REAL,DIMENSION(klon,klev) :: zrho
140	REAL, PARAMETER :: qliq=1.e-3 ! CONVECTIVE ONLY
141
142	! Je end 20140105
143
144	!JE20140724<<
145	!JE SPLA adapt:
146	!
147	!
148	!! coefficient lessivage
149	! ccntrAA_coef = 0.
150	! ccntrENV_coef = 0.
151	! coefcoli = 0.
152	!
153	!!$OMP MASTER
154	! call getin('ccntrAA_coef',ccntrAA_coef)
155	! call getin('ccntrENV_coef',ccntrENV_coef)
156	! call getin('coefcoli',coefcoli)
157	!!$OMP END MASTER
158	!!$OMP BARRIER
159	!!! JE
160	! do j=1,klev
161	! do i=1,klon
162	! imp(i,j)=0.
163	! ccntrAA_coef3d(i,j)=ccntrAA_coef
164	! ccntrENV_coef3d(i,j)=ccntrENV_coef
165	! coefcoli3d(i,j)=coefcoli
166	! enddo
167	! enddo
168	!
169	!
170	!!! for SPLA
171	!!!JEend
172	! print*,'cvltr coef lessivage convectif', ccntrAA_coef,ccntrENV_coef,coefcoli
173	!
174	!JE20140724>>
175
176	! scavON=.TRUE.
177	! if(scavON) then
178	! ccntrAA_coef = 1.
179	! ccntrENV_coef = 1.
180	! coefcoli = 1.
181	! else
182	! ccntrAA_coef = 0.
183	! ccntrENV_coef = 0.
184	! coefcoli = 0.
185	! endif
186
187	! ======================================================
188	! calcul de l'impaction
189	! ======================================================
190	!initialisation
191	do j=1,klev
192	do i=1,klon
193	imp(i,j)=0.
194	enddo
195	enddo
196	!JE init 20140103
197	!! impaction sur la surface de la colonne de la descente insaturee
198	!! On prend la moyenne des precip entre le niveau i+1 et i
199	!! I=3/4* (P(1+1)+P(i))/2 / (sigdrrho_l)
200	!! 1000kg/m3= densité de l'eau
201	!! 0.75e-3 = 3/4 /1000
202	!! Par la suite, I est tout le temps multiplié par sig_d pour avoir l'impaction sur la surface de la maille
203	!! on le néglige ici pour simplifier le code
204	! do j=1,klev-1
205	! do i=1,klon
206	! imp(i,j) = coefcoli0.75e-3/rdrop &
207	! 0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j))
208	!! rho(i,j)=pplay(i,j)/(rd*te(i,j))
209	! enddo
210	! enddo
211	!JE end 20140103
212
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).lt.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.lt.klev) then
265	if(epIN(i,j).lt.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).lt.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).lt.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).lt.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).lt.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).lt.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)).lt.1.e-10&
305	.and.Pa(i,j).lt.1.e-10.and.Pm(i,j).lt.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 .EQ. 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).gt.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).gt.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.ge.icb(i).and.j.le.inb(i)) then
530	if(clw(i,j).gt.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.ge.icb(i).and.k.le.inb(i).and.&
547	j.le.inb(i)) then
548	if(elij(i,k,j).gt.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).gt.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).gt.mp(i,j).and.mp(i,j+1).gt.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).gt.mp(i,j+1).and.mp(i,j).gt.1.e-10) then! entrainement
584	if(j.eq.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)).lt.1.e-25) then !si denominateur nul (il peut y avoir des mp!=0)
602	kappa(i,j)=1.
603	if(j.eq.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).gt.mp(i,j).and.mp(i,j+1).gt.1.e-10) then
606	qDi(i,j,it)=qDi(i,j+1,it)
607	elseif(mp(i,j).gt.mp(i,j+1).and.mp(i,j).gt.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.gt.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.eq.1) then
632	if(mp(i,2).gt.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).gt.mp(i,j).and.mp(i,j+1).gt.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).gt.mp(i,j+1).and.mp(i,j).gt.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	!------------------------------------------------------------- endif ! 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.le.inb(i).and.k.gt.1) then ! tendances dans le nuage
754	!------------------------------------------------------------- detrainement
755	if(mp(i,k+1).gt.mp(i,k).and.mp(i,k+1).gt.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).gt.mp(i,k+1).and.mp(i,k).gt.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.eq.1) then ! mp(1)=0.
789	if(mp(i,2).gt.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