Context Navigation

cvltr_spl.F90 @ 5103

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

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

Download in other formats:

Original Format